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1  Abstract 

The  objective  of  this  project  was  to  understand,  design,  and  evaluate  composite  materials 
containing  engineered  microstructures  that  display  negative  stiffness  (NS)  and  negative  Poisson's  ratio 
(auxetic)  behavior.  The  ultimate  aim  for  understanding  such  microstructural  elements  is  to  enable  the 
creation  of  composite  materials  that  are  significantly  more  dissipative  than  currently  available  materials 
(see  Figure  1)  to  enhance  existing  mechanical  absorption  and  isolation  capabilities  without  degrading 
stiffness  performance.  Progress  has  been  made  in  the  design  and  analysis  of  NS  microstructures  that 
exploit  buckled  element  geometry  to  produce  NS  behavior.  The  nonlinear  effective  material  properties 
of  the  microstructure  have  been  obtained  using  finite  element  analysis  and  the  overall  behavior  of  a 
composite  material  containing  NS  inclusions  has  been  determined  using  combined  FEA  and  analytical 
effective  medium  theory.  Auxetic  materials  work  has  produced  particulate  composite  materials 
containing  auxetic  a-crystobalite  inclusions.  The  resulting  composite  materials  have  been  characterized 
using  quasistatic,  modal  vibration,  and  ultrasonic  methods  together  with  extensive  microscopy.  The 
knowledge  gained  during  the  work  of  this  project  can  be  used  in  future  research  to  provide  a  basis  for 
the  design  of  microstructures  that  can  aid  in  absorbing  vibroacoustic  energy  to  improve  overall 
composite  material  performance  through  variations  of  engineered  microstructure. 


10* 

Metals  Proposed 


Absorption  /  Density 

Figure  1:  Plot  (log-log  scale)  of  stiffness  and  mechanical  absorption  normalized  to  mass  density  of  existing  materials  and 
composites  containing  negative  stiffness  and  auxetic  inclusions.  The  proposed  materials  lie  in  the  white  space  of  currently 
un-attainable  properties.  Estimates  were  found  using  Mori-Tanaka  model  with  2%vol  of  negative  inclusion. 
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3  International  Collaboration 


Research  supported  by  this  project  was  carried  out  as  an  international  collaboration  between  the 
University  of  Texas  at  Austin  in  the  United  States  and  the  Universities  of  Bolton  and  Bristol  in  the  United 
Kingdom.  This  project  represents  the  first  effort  to  establish  and  foster  an  international  collaboration 
between  these  centers  of  excellence  in  the  areas  of  negative  stiffness  and  negative  Poisson's  ratio 
materials.  The  principle  investigator  and  US  team  is  located  at  the  University  of  Texas  at  Austin  and  has 
led  recent  investigations  into  the  design  and  testing  of  negative  stiffness  mechanisms  for  vibration 
isolation  [la, lb].  The  members  of  the  US  team,  Drs.  Michael  Haberman,  Carolyn  Seepersad,  and  Preston 
Wilson,  have  expertise  in  the  areas  of  material  modeling,  materials  design,  wave  propagation,  and 
dynamic  material  testing  and  characterization.  The  UK  team  has  been  assembled  with  researchers  at  the 
Universities  of  Bolton  and  Bristol,  consisting  of  Drs.  Andrew  and  Kim  Alderson  at  Bolton  and  Dr.  Fabrizio 
Scarpa  at  Bristol.  These  researchers  have  extensively  studied  the  static  and  dynamic  behavior  and 
fabrication  of  auxetic  materials.  Their  combined  areas  of  expertise  include  modeling,  fabrication,  and 
characterization  of  composite  and  honeycomb  materials  for  both  static  and  dynamic  applications. 

4  Scientific  Progress  and  Accomplishments 

As  mentioned  in  the  abstract,  the  project  work  consists  of  research  in  two  coupled  areas: 
composites  containing  (/)  negative  stiffness  (NS)  and  (/#')  negative  Poisson's  ratio  (auxetic)  inclusions. 
Initial  efforts  that  incorporates  NS  elements  into  an  auxetic  network  has  also  been  investigated  during 
this  project.  Progress  on  NS  portions  of  the  research  is  presented  and  discussed  in  Section  4.1  while 
auxetic  material  research  is  discussed  in  Section  4.2.  Based  on  the  information  presented  in  those 
sections,  insight  for  associated  future  research  is  provided  in  Section  5. 

4.1  Negative  Stiffness  Materials 

The  negative  stiffness  research  effort  of  this  project  was  guided  by  the  materials  design  approach 
schematized  in  Figure  2.  The  results  presented  in  this  final  report  focus  primarily  on  steps  (/)  and  (//)  in 
the  schematized  multiscale  material  model  shown  in  the  blue  box  of  Figure  2.  The  two  steps  are 
microstructural  modeling  and  mesoscale  homogenization  of  candidate  microstructures,  respectively, 
which  are  discussed  in  Sections  4.1.2  -  4.1.6.  Progress  made  in  these  areas  has  provided  the  research 
team  with  a  comprehensive  understanding  and  robust  modeling  capability  that  captures  the  behavior  of 
candidate  microstructures  exhibiting  negative  stiffness  at  the  mesoscale.  Initial  progress  on  the 
exploration  of  the  multi-dimensional  nonlinear  design  space  using  surrogate  models  is  reported  in 
Sections  4.1.9-4.1.11.  This  tools  and  understanding  provided  by  the  coupled  multiscale  model  and 
surrogate  modeling  approach  can  be  leveraged  in  future  materials  design  work  to  tailor  microstructure 
to  produce  NS  for  improved  stiffness  and  loss  performance  at  the  macroscale.  Limited  work  has  also 
been  done  on  step  (///)  of  the  multiscale  material  model  and  is  discussed  in  Section  4.1.7. 
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Figure  2:  Schematic  of  the  multi-scale  materials  design  approach  devised  to  design  negative  stiffness  inclusions  for  improved 
macroscopic  material  performance. 

Before  discussing  the  results,  it  is  important  to  define  what  is  meant  by  the  terms  micro-,  meso-, 
and  macroscale  in  this  work  and  to  introduce  the  multiscale  modeling  scheme  at  the  core  of  the 
materials  design  effort.  In  general,  the  material  system  is  assumed  to  be  well  described  on  three  distinct 
length  scales:  the  micro-,  meso-,  and  macroscales,  in  ascending  order  and  indicated  by  boxes  (#'),  (//),  and 
(///)  in  Figure  2,  respectively.  Further,  the  materials  design  approach  assumes  that  the  overall  behavior  of 
the  composite  material,  whose  ultimate  purpose  is  to  absorb  vibro-acoustic  energy,  can  be  predicted 
using  a  multiscale  model  that  takes  structure  and  properties  at  all  length  scales  into  account.  Because 
the  microscale  behavior  required  to  produce  NS  phenomena  is  inherently  nonlinear,  the 
homogenization  scheme  must  also  admit  nonlinear  microscale  behavior  and  the  design  methodology 
must  be  of  sufficient  sophistication  to  enable  efficient,  yet  complete,  design  space  exploration  of 
nonlinear  systems.  To  achieve  these  goals,  two  distinct  scale  transition  models  were  developed  and  a 
robust  design  space  exploration  scheme  was  adopted.  The  first  scale  transition  model  is  the  micro-  to 
mesoscale  transition  model  and  the  second  is  the  meso-  to  macroscale  transition  model,  indicated  in  the 
white  boxes  labeled  (ii)  and  (iii)  of  Figure  2,  respectively.  Note  that  the  output  of  the  first  scale  transition 
model  is  part  of  the  input  to  the  second.  These  scale  transition  models  are  described  in  general  in  the 
following  paragraphs  and  in  detail  in  Sections  4. 1.5-4. 1.7.  The  design  space  exploration  is  introduced 
and  discussed  in  Sections  4.1.9-4.1.11. 

The  materials  modeled  in  this  portion  of  the  project  are  assumed  to  consist  of  a  continuous  matrix 
material  containing  a  known  volume  fraction  of  structured  inclusions.  The  macroscale  is  denoted  with 

the  length  scale,  Lu,  which  satisfies  the  condition  Lmj Lu  <C  1.  Here  Lm  is  the  descriptive  length  of 
the  structured  inclusions  and  is  the  descriptive  dimension  of  the  mesoscale.  At  the  macroscale  the  the 
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particulate  composite  material  is  assumed  to  behave  as  an  effective  homogeneous  medium.  The 
properties  of  the  effective  medium  can  be  estimated  using  well-established  effective  medium  theories 
(EMT)  if  the  elastic  properties  of  the  mesoscale  structures  as  well  as  their  volume  fraction  and  spatial 
orientation  distribution  are  known.  This  aspect  of  the  multiscale  material  model,  known  as  the  meso-to 
macroscale  transition  model,  is  addressed  in  Section  4.1.7.  In  practice  one  notes  that  a  scale  discrepancy 

of  Lmj Lu  ~  0  (l/io)  is  sufficient  to  enable  reliable  estimates  of  the  macroscopic  material  properties 

using  classical  EMT  methods  [lc]. 

In  the  present  work,  as  in  most  multiscale  modeling  literature  [lc],  the  mesoscale  is  simply  an 
intermediary  scale  between  the  inclusion  structure  feature  size  and  macroscopic  length  scale.  One  key 
focus  of  the  research  reported  here  is  the  development  of  a  nonlinear  homogenization  model  to  allow 
for  the  accurate  approximation  of  the  effective  nonlinear  elastic  response  of  an  arbitrarily  shaped 
structured  inclusion  using  finite  element  analysis  (FEA).  This  homogenization  step  is  known  as  the  micro- 
to  mesoscopic  scale  transition  model.  It  is  this  scale  transition  model,  a  very  unique  contribution  of  the 
work  associated  with  this  project,  which  that  enables  the  use  of  classical  EMT  to  predict  the  macroscopic 
properties  and  their  dependence  on  microscale  structure  and  the  associated  NS.  The  final  length  scale  of 
interest,  the  microscale,  is  defined  as  the  descriptive  lengths  of  the  structure  of  within  the  inclusions 

and  is  denoted  as  L  .  It  is  the  microstructure  that  leads  to  mesoscopic  NS  behavior  and  significant 

increases  in  absorptive  capacity  on  the  macroscopic  scale.  A  detailed  understanding  of  the  impact  of 
microscale  structure  to  produce  NS  is  therefore  the  principal  importance  for  this  project. 

4.1.1  Negative  stiffness 

Negative  stiffness  systems  are  rare  in  nature  and  therefore  not  as  well  understood  or  utilized  as 
positive  stiffness  systems.  Positive  stiffness  systems  are  simply  defined  by  the  ability  of  the  system  to 
offer  increasing  resistance  to  deformation  due  to  externally  applied  forces.  On  the  contrary,  a  negative 
stiffness  system  will  offer  decreasing  resistance  to  an  applied  load  and  can  even  assist  the  applied  force 
in  deforming  the  system  [1],  [2].  This  behavior  leads  to  significantly  larger  system  deformation  than 
would  be  observed  for  the  same  force  applied  to  a  positive  stiffness  system.  Recall  that  the  stiffness  of 
any  object,  k ,  is  the  force  required  to  deform  the  object  by  a  certain  distance,  or  F  =  kAx  .  One  very 
simple  example  of  a  system  displaying  negative  stiffness  properties  is  a  transversely  loaded  buckled 
beam.  A  buckled  beam  is  a  simple  example  of  a  bistable  system.  Negative  stiffness  can  be  observed  for 
such  a  bistable  system  in  one  of  two  scenarios,  the  first  of  which  is  depicted  by  Figure  3.  In  this  case,  the 
beam  starts  in  one  of  the  two  stable  configurations  and  driven  to  the  other  bistable  configuration.  In  this 
scenario,  the  buckled  beam  initially  resists  deformation,  but  once  a  certain  force  threshold  is  reached, 
the  beam  "snaps  through"  to  a  new  configuration  and  remains  in  that  configuration  even  after  the  force 
is  removed.  In  this  case,  the  negative  stiffness  behavior  occurs  during  the  snap  through,  is  transient,  and 
leads  to  large,  and  permanent,  deformations.  The  second  scenario  where  negative  stiffness  can  be 
observed  is  the  case  where  a  buckled  beam  is  constrained  at  its  center  point  by  another  element,  such 
as  a  spring.  In  this  case,  the  beam  remains  in  position  (2)  of  Figure  3  but  will  move  when  perturbed  by  an 
outside  force.  The  constraining  element  will  then  force  the  beam  to  return  to  position  (2)  when  the 
outside  force  is  removed.  This  scenario  represents  the  case  of  constrained  negative  stiffness  (CNS), 
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which  is  very  useful  for  vibroacoustic  damping  where  loads  are  low  enough  in  amplitude  to  not  force  the 
beam  to  oscillate  between  bistable  positions. 


(D  (2)  (3) 


F(Xj) 


Figure  3:  Schematic  of  the  behavior  of  a  beam  of  orignial  length  L0  with  buckled  length  L0  -  6 L.  The  beam  is  loaded  in 
displacement  control  from  one  bistable  position,  (configuration  1),  to  the  other,  (configuration  3). 

The  behavior  of  both  of  these  scenarios  and  how  it  relates  to  negative  stiffness  can  be  related  to  the 
positions  and  forces  shown  in  Figure  3.  Notably,  the  buckled  beam  force  versus  displacement  response 
is  in  stark  contrast  to  that  of  an  unbuckled  beam  where  an  applied  transverse  force  will  cause  the  beam 
to  deflect  in  proportion  to  the  applied  external  force  only  to  return  to  its  original  position  when  all 
forces  are  removed. 

To  better  understand  how  this  behavior  is  indicative  of  negative  stiffness  and  what  the  roots  of  this 
behavior  are  in  terms  of  mechanical  metrics,  it  is  useful  to  consider  force  and  internal  strain  energy 
versus  displacement  behavior.  Both  of  these  metrics  are  depicted  in  curves  shown  in  Figure  4  and 
compared  with  the  behavior  of  an  unbuckled  beam  for  reference.  The  left  panel  of  this  figure  shows  the 
normalized  force  versus  displacement  of  the  beam  center  point.  The  normalized  force  is  given  as 
2 F  /  Pcr  where  F  is  the  transversely  applied  load  and  Pcr  is  the  axial  load  that  induces  buckling  as 
predicted  by  Euler-Bernouilli  beam  theory.  The  normalized  displacement  is  Ax  /  L  where  Ax  is  the 

displacement  from  the  x  =  0  position  and  L  is  the  free  beam  length.  There  are  several  points 
illustrated  in  these  curves  that  are  worth  highlighting.  First  of  all,  note  that  positions  (1),  (2),  and  (3)  all 
represent  configurations  where  zero  force  is  required  to  maintain  that  shape.  In  other  words,  theory 
predicts  that  the  c-shape  of  configuration  (points  1  and  3)  and  the  s-shape  configuration  (point  2)  are 
points  of  'stability.'  These  stability  conditions  can  be  qualified  through  inspection  of  the  internal  energy, 
U(x),  versus  displacement  curve  shown  in  the  right  panel.  Namely,  these  stability  positions  correspond 
to  points  of  zero  slope  on  the  strain  energy  curve,  i.e.  dUjdx  =  0.  This  is  not  surprising  since 

mechanics  relates  internal  energy,  force,  F(x),  and  stiffness,  k(x),  through  the  relations  [3] 


and 


dU(x) 

dx 


(1) 
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(2) 


d2u(x) 

dx2 


-0.01  -0.005  0  0.005  0.01 


AX/2L0 

(a) 


Figure  4:  (a)  Force  versus  displacement  and  (b)  internal  (strain)  energy  for  a  transversely  loaded  buckled  beam  such  as  the 
one  shown  in  Figure  3. 


Though  these  three  configurations  can  be  qualified  as  points  of  stability,  physical  intuition  suggests  that 
if  one  were  to  place  the  beam  in  position  (2)  it  would  not  remain  there  when  constraining  forces  are 
removed.  This  is  a  result  of  the  fact  that  the  position  is  not  unconditionally  stable,  specifically,  it  is  in  a 
position  of  unstable  equilibrium.  For  stable  equilibrium,  the  change  in  the  internal  energy  resulting  from 
perturbations  about  the  equilibrium  position  must  be  positive.  Mathematically  this  means  that  the 

curvature,  and  therefore  the  stiffness,  of  the  internal  energy  curve  must  be  positive,  d2U J dx2  >  0. 

This  requirement  implies  that  any  deformation  from  the  position  leads  to  an  increase  in  internal  energy, 
which  is  another  way  of  stating  that  a  stable  equilibrium  position  is  one  where  the  strain  energy  is  at  a 
local  minimum,  also  known  as  an  energy  well  [4].  Using  this  criteria,  unconditional  stability  occurs  at 

positions  (1)  and  (3)  where  the  slope  of  the  internal  energy  curve  is  zero,  dUjdx  =  F^xj  =  0,  and 

where  its  curvature  is  positive,  d2U ^ j  dx 2  =  k  (a;)  >  0 .  The  unstable  equilibrium  position,  configuration 

(2),  is  characterized  as  a  point  of  zero  force  where  the  curvature  of  the  internal  energy  curve  is  negative, 
i.e.  a  local  maximum  in  the  strain  energy  of  the  beam.  Any  unconstrained  beam  in  this  configuration  will 
immediately  be  driven  to  one  of  the  two  unconditionally  stable  positions  by  any  minor  perturbation. 
Another  very  interesting  point  shown  by  Figure  4  is  that  the  force-displacement  curve  for  the  buckled 
beam  is  non-monotonic.  This  implies  that  the  slope,  and  therefore  the  stiffness,  of  the  transverse  force 
versus  displacement  of  the  buckled  beam  will  display  negative  stiffness.  In  other  words,  the  existence  of 
a  bistability  alone  implies  that  negative  stiffness  behavior  exists  and  may  potentially  be  advantageously 
used. 

The  force  and  internal  energy  versus  displacement  curves  very  clearly  illustrate  why  the  two 
scenarios  discussed  above  will  display  negative  stiffness.  In  the  first  case,  when  the  beam  is  driven  from 
one  bistability  to  the  other,  the  stiffness  offered  by  the  beam  is  negative  between  inflection  points  of 
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the  internal  energy  curve.  This  region  is  highlighted  in  the  left  panel  of  Figure  4.  In  the  second  scenario, 
the  beam  is  constrained  to  position  (2),  which  is  a  local  maximum  in  strain  energy.  As  was  pointed  out 
above,  a  local  maximum  in  internal  energy  corresponds  to  a  region  of  negative  curvature  in  the  energy 
curve.  Because  the  stiffness  of  the  element  is  equal  to  the  curvature  of  the  internal  energy  curve,  the 
beam  will  categorically  display  negative  stiffness  in  the  region  for  all  perturbations  about  its  constrained 
position.  The  constrained  behavior  of  the  system,  referred  to  above  as  CNS,  is  not  limited  to  large 
displacements  and  is  therefore  of  high  interest  for  vibroacoustic  damping  applications. 

It  is  the  focus  of  this  portion  of  the  project  to  exploit  the  unique  non-monotonic  force-displacement 
nature  of  bistable  systems  to  produce  highly  absorptive  composite  materials.  This  can  be  achieved  by 
designing  and  then  embedding  negative  stiffness  inclusions  in  a  continuous  viscoelastic  matrix  to 
effectively  enhance  the  damping  properties  of  the  matrix  material.  As  such,  the  objectives  of  this 
portion  of  the  research  were  to  (/')  discover  microstructure(s)  that  displays  negative  stiffness  and  (/'/) 
validate  a  design  approach  using  finite  element  analysis  (FEA)  that  has  the  ability  to  quantify  various 
microstructures  that  yield  negative  stiffness  and  approximate  their  behavior  as  an  effective  stiffness 
tensor. 

4.1.2  Mesoscale  inclusion  characterization 

Though  modeling  conceptually  begins  at  the  microscale,  it  is  practically  useful  to  begin  with  a 
discussion  of  the  mesoscale  inclusion  properties.  Namely,  it  is  of  interest  to  detail  how  to  reliably 

approximate  a  mesoscopic  effective  stiffness  tensor,  C'"jd ,  of  an  inclusion  with  arbitrary  internal 

geometry  based  on  its  stress-strain  (force-displacement)  behavior  when  subjected  to  various  loading 
conditions.  With  a  reliable  method  to  approximate  mesoscopic  stiffness  in  place,  it  is  then  possible  to 
focus  on  direct  analysis  and  design  of  the  inclusion  microstructure.  The  work  carried  out  during  this 
project  employed  two  approaches  to  calculate  the  effective  mesoscopic  stiffness.  The  first  is  referred  to 
as  the  direct  energy  method  which  is  reliable  for  linear  elastic  materials  and  structures.  The  second  is 
called  the  energy  derivative  method  which  is  well  suited  for  materials  and  structures  that  display 
nonlinear  constitutive  stress-strain  relationships. 

The  method  referred  to  as  the  direct  energy  method  is  one  that  was  employed  by  Odegard  to 
calculate  the  effective  behavior  of  a  representative  volume  element  of  multiphase  continuous 
piezoelectric  materials  using  FEA  [5].  The  methodology  is  straightforward  and  involves  computing  the 
change  in  internal  energy  within  a  material  volume,  known  as  the  representative  volume  element  (RVE), 
which  is  assumed  to  represent  the  overall  response  of  the  material  of  interest  when  the  volume  is 
subjected  to  a  prescribed  set  of  boundary  displacements.  The  mathematical  details  are  summarized  by 
Eqs.  (3)-(7)  and  Table  1.  Using  the  Einstein's  summation  convention  for  indicial  notation,  the  linear 

elastic  constitutive  equation  relating  stress,  a..,  and  strain,  Eu,  in  stiffness  form  is 


a  — 

ij  ijld  kl 


(3) 
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where  C.jkl  is  the  stiffness  tensor  of  the  material  (tacitly  assumed  to  be  linear  in  this  formulation)  and 
the  (Green's)  strain  at  a  material  point  defined  as  the  symmetric  part  of  the  gradient  of  the 
displacement  field,  Eu  =  (ukl  +  ulk)/ 2  •  The  change  in  specific  internal  energy  at  a  material  point, 

U® ,  due  to  an  imposed  strain  is  given  by 

u? 4<wv  <4) 

Note  that  the  middle  equality  in  Eq.  (4)  is  valid  for  materials  with  arbitrary  elastic  stress-strain  relations 
while  the  approximate  equality  on  the  far  right-hand  side  (RHS)  is  true  for  linear  elastic  materials  for  all 
imposed  strains  or  nonlinear  elastic  materials  in  the  limit  where  the  imposed  displacements  at  the 
boundaries  are  small.  For  orthotropic  materials,  Eq.  (3)  can  be  represented  as  follows  with  Voigt 
notation 
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where  E .  £..  indicating  small  strains  and  7  =  2 e..  to  preserve  symmetry. 

IJ  IJ  1  IJ  l] 

Using  COMSOL  Multiphysics  finite  element  modeling  software,  homogeneous  strains  are  applied  at 
boundaries  of  a  cube  of  material.  At  boundary,  B ,  the  displacements,  u.  (B^J,  can  be  related  to  the 

imposed  strain,  e° ,  and  the  dimension  of  the  RVE  cube,  L  ,  in  a  general  way  by 

ij  m 


£  L  n., 

ijmj 


(6) 


where  n  is  the  unit  normal  vector  of  the  RVE  face  where  the  displacement  is  imposed.  Note  that  this 
formulation  has  made  the  assumption  that  the  strain  in  the  RVE  is  well  represented  as  being 
homogeneous  within  the  RVE.  The  resulting  total  elastic  strain  energy  is  then  defined  by  the  following 
relation: 


[7  =  YK  =  —C.ve0.el, 

E  /  j  E  2  M 

m= 1 


(7) 
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where  U™  is  the  elastic  strain  energy  of  element  m,  n  is  the  total  number  of  finite  elements  and  V  is  the 

volume  of  the  cube.  The  total  energy  is  then  summed  over  the  entire  volume.  In  order  to  determine  the 
elements  of  the  stiffness  tensor  it  suffices  to  work  out  the  mathematical  relations  of  specific  applied 
boundary  conditions  and  the  resulting  strain  energies  by  dividing  the  strain  energy  by  the  imposed  strain 
and  a  constant.  Table  1  summarizes  the  boundary  conditions  and  related  strain  energy  relations  for  the 
nine  independent  constants  of  an  orthotropic  material. 


Table  1:  Summary  of  the  displacements  and  expressions  to  determine  the  independent  constants  of  an  orthotropic  medium 
using  the  direct  energy  approach. _ 
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Property 

Applied  strain 

Displacements 

Elastic  energy 
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Research  during  this  project  has  illustrated  that  the  direct  energy  method  is  ill-suited  to  determine 
the  stiffness  of  structures  displaying  negative  stiffness.  There  are  two  primary  reasons  for  this 
shortcoming.  First,  the  direct  energy  approach  is  formulated  such  that  minor  errors  in  the  calculated 
strain  energy  can  result  in  very  large  errors  in  the  stiffness  coefficient  estimate  because  of  the 
requirement  of  small  imposed  strains  and  the  presence  of  the  strain  in  the  denominator  of  the  stiffness 
calculation.  Second,  if  the  material  or  microstructure  exhibits  nonlinear  stress-strain  behavior,  this 
approach  has  the  potential  to  produce  non-physical  stiffness  results.  Given  that  the  focus  of  the  current 
research  is  the  determination  of  inclusion  structures  that  display  negative  stiffness,  and  that  negative 
stiffness  requires  the  presence  of  bistability,  which  is  inherently  nonlinear,  the  direct  energy  approach  is 
likely  to  provide  erroneous  estimates  of  the  stiffness  of  the  structures  of  interest. 

Fortunately,  a  method  that  uses  much  the  same  approach,  what  developed  by  the  PI  and  has  been 
named  the  energy  derivative  method.  This  approach  can  be  implemented  to  estimate  the  effective 
nonlinear  behavior  of  a  structure  that  undergoes  finite  deformation  and  does  not  have  the  weaknesses 
of  the  direct  approach.  Much  like  the  direct  approach  described  above,  the  energy  derivative  method 
also  applies  prescribed  strains  on  the  surface  of  the  RVE  and  calculates  the  resulting  strain  energy.  The 
primary  difference  is  that  the  energy  derivative  method  requires  the  calculation  of  the  strain  energy  for 
a  large  range  of  imposed  deformations.  Once  the  resulting  strain  energies  are  calculated,  the  stiffness  of 
the  element  at  a  specified  strain  can  be  calculated  from  the  determination  of  the  local  curvature  of  the 
strain  energy  versus  strain  plot.  This  approach  does  not  require  the  stress-strain  relation  to  be  linear 
thereby  implying  that  only  the  very  general  center  equality  of  Eq.  (4)  hold.  Research  conducted  during 
this  project  has  shown  that  the  energy  derivative  method  is  a  robust  approach  for  determining  local 
stiffness  and  eliminates  the  error  amplification  inherent  in  the  direct  approach.  Equations  (8)-(10) 
describe  the  relationships  between  the  strain  energy  curves  created  by  applying  a  range  of  deformations 
at  the  boundaries  of  the  RVE  similar  to  those  given  in  Table  1  and  the  nonlinear  stiffness  curves  for  the 
mesoscopic  structure.  Specifically,  the  total  strain  energy  is  related  to  the  assumed  homogeneous  stress 
and  strain  fields  as 


uE=iy; is  ie 

m—l 


(8) 
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This  implies  that  the  stress  at  any  strain  level  can  be  computed  from  the  derivative  of  the  strain  energy  - 
strain  curve: 


a  . 

y 


V  de 

y 


(9) 


Finally,  the  local  stiffness  is  then  simply  found  as  the  curvature  of  the  strain  energy  versus  strain  curve, 


,  .  da..(E.. 

C.v  ” 

y ki  \  y/  /)p 

utz'ki 


V  dE 2 


(10) 


Equation  (10)  implies  that  the  prescribed  displacements  in  Table  1  can  be  applied  to  an  arbitrary  RVE  for 
varying  magnitudes  of  displacement  and  the  resulting  total  elastic  strain  energy  in  the  RVE  can  be 
computed  due  to  these  imposed  strains.  Calculations  of  the  local  curvature  at  any  given  displacement 
(or  strain)  can  then  be  related  to  the  corresponding  inclusion  stiffness  through  Eq.  (10).  This  approach  is 
first  benchmarked  in  Section  4.1.3  and  then  employed  to  calculate  the  effective  mesoscopic  stiffness  of 
a  candidate  microstructure  in  Section  4.1.4. 

4.1.3  Benchmarking  the  energy  derivative  approach 

To  benchmark  the  energy  derivative  approach,  FEA  was  carried  out  on  a  solid  block  of  isotropic 
steel  (E  =  200  GPa,  v  =  0.33).  The  following  figures  contain  the  respective  strain  energy-displacement, 

force-displacement  and  stiffness-displacement  values  for  Cn,  C12  and  Cm .  Neglecting  numerical 

round-off  from  derivative  calculations,  the  values  converge  to  296  GPa,  146  GPa  and  75.2  GPa, 
respectively,  for  the  full  range  of  imposed  displacement.  These  values  are  equivalent  to  the  known 
properties  of  isotropic  steel,  providing  a  good  initial  benchmark  of  the  energy  derivative  approach. 
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Figure  5:  Calculated  strain  energy,  force,  and  stiffness  as  a  function  of  displacement  imposed  on  the  RVE  for  a  steel  cube 
using  the  energy  derivative  approach.  For  steel  Cn  =  A  +  2p  =  296  GPa,  C12  =  A  =  146  GPa  and  C66  =  p  =  75.2  GPa. 

To  provide  a  more  robust  benchmark  of  the  energy  derivative  approach,  it  is  of  interest  to  model 
the  behavior  of  a  heterogeneous  system  and  compare  the  FEA  results  to  results  from  well-established 
EMT.  This  was  accomplished  by  modeling  a  RVE  consisting  of  a  cube  of  steel  containing  a  spherical  void 
of  varying  volume  fraction.  The  same  process  for  obtaining  stiffness  values  was  performed  on  this  RVE 
as  for  the  homogeneous  cube  and  the  results  are  displayed  in  Figure  6.  This  figure  contains  three  plots, 
one  of  each  constant  of  interest  (Cn,  C12  and  Cm )  for  the  heterogeneous  RVE  as  the  void  fraction  is 

varied  from  0  to  20%  by  volume.  The  values  are  normalized  by  the  element  of  the  stiffness  tensor  of 
steel.  The  curves  clearly  show  very  strong  agreement  with  the  analytical  models  which  provides  a  very 
reliable  benchmark  of  the  energy  derivative  approach  for  approximating  the  overall  properties  of  a 
heterogeneous  RVE.  With  a  method  for  approximating  the  effective  properties  of  a  structured  inclusion 
in  place,  attention  is  now  turned  to  modeling  microstructural  elements  that  display  negative  stiffness. 
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Figure  6:  Effective  properties  of  a  RVE  of  steel  containing  varying  amounts  of  spherical  void  (see  inset  in  top  panel).  Stiffness 
values  were  calculated  in  FEA  using  the  derivative  approach  and  are  compared  with  analytical  estimates  using  the  Self- 
Consistent  (SC),  Differential  Effective  Medium  (DEM),  and  Mori-Tanaka  (MT)  models  along  with  the  Hashin-Shtrikman  Upper 
(HS+)  and  lower  (HS  )  bounds. 

4.1.4  Micro  structure  modeling 

As  mentioned  at  the  outset  of  this  section,  the  physics  underlying  negative  stiffness  behavior  that  is 
of  interest  to  this  work  is  found  and  easily  illustrated  through  an  investigation  of  a  buckled  beam.  For 
this  reason,  it  is  of  interest  to  investigate  a  buckled  beam  through  both  analytical  and  FE  methods.  This 
serves  the  dual  purpose  of  providing  a  clear  example  of  the  desired  system  behavior  and  a  benchmark  of 
FE  modeling  of  post-buckled  beam  behavior.  Consider  the  system  shown  in  Figure  7  consisting  of  a  beam 
with  fixed  ends  subjected  to  both  axial  and  transverse  forces  labeled  as  N  and  P,  respectively. 


I 

Figure  7:  Geometry  and  variables  for  analysis  of  the  behavior  of  pre-  and  post-buckled  beam. 

Through  bending  analysis,  it  can  be  shown  that  the  transverse  beam  stiffness  at  the  center  of  the  beam, 
kb,  is  equal  to: 
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Pk 

tan^kL/^-kL/i 


(11) 


where  k  =  j El  ,  E  is  the  Young's  modulus  of  the  beam  materials,  and  /  is  the  bending  moment  of 

inertia  of  the  beam.  One  can  also  readily  plot  transverse  stiffness  versus  axial  compression,  shown  in 
Figure  8,  by  relating  the  axial  compression  force,  N,  to  the  resulting  compressive  displacement  along  the 
beam  axis,  u,  with  u  =  NLj  AE . 

Equation  (11)  and  the  relations  above  can  be  used  to  validate  FEA  results  for  a  2D  plane-strain 
beam.  A  simple  beam  model  was  created  for  an  alumina  (f  =  386  GPa,  v  =  0.22)  beam  of  length  .0017m 
and  height  of  50  pm.  A  point  displacement  of  1  pm  was  applied  at  the  center-top  of  the  beam  at  the 
same  time  as  the  beam  ends  were  subjected  to  varying  axial  compressive  displacements.  The  reaction 
force  in  the  y-direction  was  monitored  at  the  point  of  imposed  displacement  and  the  transverse  stiffness 

of  the  system  at  the  center  was  then  calculated  by  noting  that  kEEA  m  F™eas  j  ul™ posed  .  Transverse 

stiffness  values  obtained  using  FEA  and  Eq.  (11)  are  plotted  in  Figure  8  for  varying  amounts  of  axial 
compression.  The  results  indicate  a  very  good  agreement  between  the  analytical  and  FEA  models.  The 
slight  deviations  are  attributed  to  discrepancies  between  the  analytical  equation  which  assumes  a  ID 
beam  and  the  2D  plane-strain  approximation  used  in  FEA.  For  the  purposes  of  this  research,  these 
results  validate  this  FEA  approach  for  modeling  pre-  and  post-buckled  beam  behavior.  It  is  noteworthy 
that  the  stiffness  curve  crosses  the  point  of  zero  stiffness  for  an  axial  compression  of  480  nm.  This  value 
is  equivalent  to  the  critical  compression  load,  iV  ,  i.e.  the  loading  at  which  buckling  occurs,  and  is 
relevant  to  the  discussion  of  the  results  that  are  presented  below. 
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To  further  explore  the  FEA  model,  several  studies  were  performed  in  order  to  validate  the  ability  of 
the  FE  model  to  capture  negative  stiffness  behavior.  Specifically,  the  beam  was  subjected  to  3  levels  of 

compression:  u  —  0,  u  —  u  >  ?/critical  and  the  results  are  shown  in  Figure  9.  In  the  first  case, 

the  beam  acts  as  a  conventional  positive  stiffness  system.  That  is,  the  transverse  force  versus 
displacement  curve  is  monotonically  increasing  and  the  strain  energy  has  a  minimum  value  at  Ax  =  0. 
This  is  reflected  in  the  L  —  LQ  case  shown  in  Figure  9.  As  the  amount  of  compression  increases  to  the 

point  of  criticality,  the  stiffness  approaches  zero  for  Ax  —  0.  This  is  shown  in  Figure  9  as  the  slope  of 
the  force-displacement  curve  for  the  N  =  iVcritical  case  is  shown  to  be  nearly  zero  for  a  range  of 

displacement  values  around  Ax  =  0  and  the  strain  energy  has  a  broad  range  of  displacements  around 
Ax  =  Ofor  which  it  is  nearly  constant.  A  buckled  beam  behaves  quite  different  than  both  of  the 
previous  configurations.  Starting  at  the  first  zero  crossing  on  negative  displacement  part  of  the  curve  of 
Figure  9,  the  beam  initially  resists  a  transverse  displacement  by  requiring  increasing  force  to  impose  a 
displacement.  As  the  displacement  increases  to  a  critical  value,  the  beam  will  show  decreasing 
incremental  changes  in  force  for  an  increase  in  displacement  until  it  reaches  the  second  zero  crossing  in 
the  positive  range  of  the  imposed  displacement.  This  behavior  is  indicated  by  the  double-well  nature  of 
the  strain  energy  curve  for  the  'Buckled'  case  in  Figure  9.  The  results  for  these  three  different  axial 
loading  conditions  are  in  very  good  agreement  with  the  discussion  about  negative  stiffness  systems  in 
Section  4.1.1  and  the  FEA  approach  is  therefore  validated  for  finite  deformation  systems  that  display 
negative  stiffness. 
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Figure  9:  FEA  results  of  the  Force  (top)  and  strain  energy  per  beam  width  (bottom)  versus  displacement  curves  for  a  beam 
with  varying  axial  load. 

As  mentioned  in  Section  4.1.1,  the  fact  that  negative  stiffness  implies  the  presence  of  a  local 
increase  in  strain  energy,  they  can  only  exist  if  constrained  within  an  overall  positive  stiffness  system 
which  allows  energy  to  be  'stored'  locally  in  the  negative  stiffness  system.  When  perturbed  from  this 
constrained  position,  NS  systems  can  release  more  energy  into  their  constraining  systems  than  the 
energy  input  by  external  sources.  This  is  the  key  point  of  the  use  of  NS  elements  to  absorb  energy: 
negative  stiffness  systems  can  be  used  to  perform  work  on  their  external  constraining  system.  If  the 
surrounding  medium  is  lossy,  then  the  presence  of  a  NS  domain  will  provide  more  energy  to  work  that 
surrounding  medium  and  thereby  increase  the  absorptive  capacity  of  the  overall  system.  This  is  the  basis 
of  hypotheses  of  extreme  damping  using  negative  stiffness  systems  [6]. 

4.1.5  Mesoscale  inclusion  modeling  and  design 

The  current  inclusion  design  mimics  a  buckled  beam  system  on  four  of  the  six  faces  of  a  cube  to 
induce  negative  stiffness  in  two  orthogonal  directions.  The  structure  shown  in  Figure  10  was  designed  to 
take  advantage  of  a  fabrication  process  known  as  micro  co-extrusion  which  could  be  leveraged  to 
fabricate  and  employ  a  coefficient  of  thermal  expansion  mismatch  between  Al203  (alumina)  and  Y203 
Partially  Stabilized  Zirconia  (YTZP),  two  ceramics  that  can  be  polymerized  for  co-extrusion,  to  produce 
buckled  elements.  Given  a  sufficient  change  in  temperature,  ample  axial  compression  can  be  generated 
to  buckle  the  alumina  'web'  elements  depicted  in  Figure  10.  Attached  to  the  web  is  a  YTZP  interface  to 
the  matrix  material  which  will  surround  the  mesoscale  inclusions.  The  YTZP  interface  acts  as  a 
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distributed  to  point  force  translator  to  ensure  that  as  much  external  loading  is  applied  directly  to  the 
center  of  the  web  as  possible.  Without  this  concentration  of  force  at  the  center  of  the  web,  FEA  has 
shown  that  the  web  can  contort  into  secondary  modes,  limiting  or  even  eliminating  its  NS  response. 


SOpm 


2mm 


2mm 

Figure  10:  Mesoscale  inclusion  with  microscale  structure  that  leads  to  negative  stiffness  in  two  dimensions.  The  inclusion  is 
assumed  to  be  fabricated  from  YTZP  and  alumina  materials  that  employ  thermal  expansion  mismatches  to  induce  buckling. 

The  objective  for  the  design  of  the  microstructure  is  ultimately  to  tailor  design  parameters  for  an 
existing  application  while  considering  existing  manufacturing  processes  for  future  implementation.  This 
means  that  the  trade-space  will  be  constrained  not  only  with  regard  to  material  selection  (for  instance, 
materials  similar  to  the  ceramics  indicated  above),  but  geometry  as  well.  Because  50  pm  is  considered 
the  smallest  dimension  that  can  be  readily  manufactured  using  micro  co-extrusion,  the  initial  design  has 
restricted  the  smallest  feature  size  to  be  a  maximum  of  50  pm.  As  of  the  termination  of  this  project,  no 
further  manufacturing  processes  have  been  considered.  However,  it  is  worth  noting  that  a  robust 
multiscale  model  coupled  with  the  design  space  exploration  methods  described  in  Section  4.1.11  permit 
optimization  of  negative  stiffness  behavior  through  the  alteration  of  inclusion  structure  rather  than 
focusing  on  a  specific  manufacturing  technique.  This  opens  the  door  to  evaluating  what  can  be  achieve 
within  the  limits  of  a  specific  manufacturing  process  rather  than  searching  for  a  design  process  that  can 
achieve  the  performance  predicted  by  a  specific  geometric  configuration. 

4.1.6  FEA  of  candiate  negative  stiffness  mesoscale  inclusion 

Characterization  of  the  temperature  depending  effective  elastic  constants  of  structured  inclusions 
like  the  one  shown  in  Figure  10  requires  the  implementation  of  a  multi-step  FEA  study  that  considers 
both  the  thermal  and  mechanical  deformation.  The  process  devised  requires  two  discrete  steps.  The 
first  step  is  to  solve  the  free  thermal  deformation  for  any  given  temperature  change.  The  resulting 
deformed  state  of  the  current  inclusion  structure  is  shown  in  the  top  left  image  of  Figure  11.  The 
resulting  displacements  are  monitored  at  the  external  boundaries  of  the  inclusion.  These  displacements 
signify  the  new  thermal  equilibrium  positions  of  the  inclusion  boundaries.  The  second  step  then 
simultaneously  solves  for  the  thermal  deformation  and  the  prescribed  strain  as  outlined  in  Table  1  in 
order  to  determine  each  independent  elastic  constant  of  the  mesoscale  inclusion.  The  symmetries 
contained  in  the  geometry  lead  to  the  observation  that  Cn  =  C22,  C  =  C23 ,  and  C44  =  C55  leaving 
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the  only  independent  constants  to  determine  as  Cn,  C33,  C12,  C23 ,  Cu,  and  C&, .  Figure  12  -  Figure 

18  show  the  results  of  this  analysis  for  all  independent  constants  of  the  inclusion  while  Figure  10  shows 
the  deformed  fields  for  a  few  representative  loading  conditions  to  aid  in  understanding  the  process. 


Plane  Strain  Thermal  Contraction 


Plane  Strain  Bulk  Modulus 


Figure  11:  Representative  deformed  microstructures  showing  the  results  of  the  thermal  contraction-induced  buckling  and 
imposed  displacements  to  determine  strain  energy  curves  for  effective  mesoscopic  constant  calculation.  The  deformed 
geometric  is  greatly  exagerated  for  easy  visualization. 


The  first  set  of  results,  shown  in  Figure  12,  depicts  the  normalized  strain  energy  and  effective 
sitffness  resulting  from  a  plane-strain  area  change  (equal  displacements  along  orthogonal  directions)  for 
differing  imposed  temperature  changes.  The  stiffness  value  resulting  from  this  set  of  boundary 
conditions  is  the  2D  bulk  modulus,  K  D .  The  strain  energy  curves  for  large  temperature  change  show 

two  energy  wells,  clearly  incidcating  bistable  behavior.  Consequently,  the  lower  panel  shows  that  for  a 
range  of  imposed  displacements  the  effective  2D  bulk  modulus  will  be  negative.  It  is  worth  noting  that 
the  strain  energy  curve  is  not  symmetric  about  the  local  maximum.  This  is  due  to  the  fact  that  the 
inclusion  geometry  is  such  that  the  'buckled-out'  configuration  (shown  in  Figure  11)  has  lower  energy 
levels  than  does  the  'buckled-in'  configuration.  This  does  not  change  the  fact  that  NS  will  be  observed 
as  is  evidenced  by  the  lower  panel  in  Figure  12. 
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Figure  12:  Effective  2D  plane  strain  bulk  modulus  of  the  candiate  microstructure  for  various  imposed  displacements  for 
various  changes  in  temperature.  Top  panel  shows  the  strain  energy  curves  (normalized  by  their  maximum  value)  and  the 
bottom  panel  shows  the  bulk  modulus. 

The  results  shown  in  Figure  12  suggest  that  the  current  inclusion  design  will  indeed  show  negative 
stiffness  given  a  certain  level  of  temperature  change  to  induce  buckling  in  the  'web'  elements.  This  is 
verified  with  the  plots  of  the  remaining  effective  stiffnesses  in  Figure  13  -  Figure  18.  Specifically,  the 

constants  Cn  =  C22  and  C12  display  negative  stiffness  behavior  for  a  relatively  large  range  of  imposed 
displacements. 
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Figure  13:  Effective  plane  wave  modulus  of  the  candidate  microstructure  along  the  xx  (and  x2)  direction  at  various  imposed 
displacements  and  changes  in  temperature.  Top  panel  shows  the  strain  energy  curves  (normalized  by  their  maximum  value) 
and  the  bottom  panel  shows  the  modulus. 
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Figure  14:  Effective  plane  wave  modulus  of  the  candidate  microstructure  along  the  x2  direction  at  various  imposed 
displacements  and  changes  in  temperature.  Top  panel  shows  the  strain  energy  curves  (normalized  by  their  maximum  value) 
and  the  bottom  panel  shows  the  modulus. 
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Figure  15:  Effective  dilatation  modulus  of  the  candidate  microstructure  in  Xi-x2  plane,  C12,  at  various  imposed  displacements 
and  changes  in  temperature.  As  is  shown  in  Table  1,  this  modulus  is  determined  from  the  results  of  other  moduli,  so  the 
strain  energy  curve  is  not  shown. 
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Figure  16:  Effective  dilatation  modulus  of  the  candidate  microstructure  in  Xi-x3  plane,  C13,  at  various  imposed  displacements 
and  changes  in  temperature.  As  is  shown  in  Table  1,  this  modulus  is  determined  from  the  results  of  other  moduli,  so  the 
strain  energy  curve  is  not  shown. 
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Figure  17:  Effective  shear  modulus  of  the  candidate  microstructure  in  x2-x3  plane, ,  C44,  at  various  imposed  displacements 
and  changes  in  temperature.  Top  panel  shows  the  strain  energy  curves  (normalized  by  their  maximum  value)  and  the  bottom 
panel  shows  the  modulus. 
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Figure  18:  Effective  shear  modulus  of  the  candidate  microstructure  in  Xi-x2  plane, ,  C66,  at  various  imposed  displacements 
and  changes  in  temperature.  Top  panel  shows  the  strain  energy  curves  (normalized  by  their  maximum  value)  and  the  bottom 
panel  shows  the  modulus. 

The  above  results  validate  the  approach  of  exploiting  a  coefficient  of  thermal  expansion  mismatch 
to  obtain  a  negative  stiffness  inclusion.  It  is  important  to  note  that  the  current  level  of  temperature 
change  required  is  too  large  for  practical  application,  however  we  have  established  a  benchmarked  and 
validated  approach  for  obtaining  temperature-dependent  inclusion  stiffness  values.  Future  work 
includes  performing  optimization  of  the  inclusion  structure  and  materials  in  order  to  bring  the  required 
temperature  change  down  to  a  level  that  is  practical  for  implementation  using  known  processing 
techniques.  Structure  that  displays  negative  stiffness  but  does  not  rely  on  temperature  change,  studied 
in  micro-switch  research,  will  also  be  investigated  to  remove  that  additional  manufacturing  difficulty. 

4.1.7  Mesoscale  me  so  scale  effective  medium  modeling 

The  remaining  modeling  step  to  complete  the  multi-scale  model  of  our  structured  lossy  composite 
material  is  the  meso  maco  homogenization  modeling  step.  This  is  completed  by  implementing  an 
anisotropic  EMT  analysis  that  uses  as  an  input  the  effective  stiffness  tensor  of  the  mesoscale  inclusion 
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introduced  in  the  Sections  4.1.5  and  4.1.6.  The  micromechanical  model  employed  is  a  three-phase 
coated  inclusion  self-consistent  (SC)  EMT  developed  in  previous  work  [7].  The  model  employs 
localization  and  homogenization  of  the  effects  of  viscoelastic  contrast  due  to  the  presence  of  microscale 
coated  inclusions,  as  indicated  by  the  schematic  in  Figure  19.  The  model  can  consider  inclusions  of 
arbitrary  anisotropic  viscoelastic  behavior,  geometry,  and  orientation  distribution  to  get  a  good 
approximation  of  the  overall  viscoelastic  response  of  the  medium  and  has  been  shown  to  have  strong 
agreement  with  experimental  data  in  static  and  quasistatic  (inclusion  size  is  much  larger  than 
propagating  wavelengths)  domains  [7],  [8].  The  effective  (complex  valued)  stiffness  tensor  of  a 
composite  containing  coated  inclusions  are  computed  using  Eqs.  (12)  -  (16). 


RVE  Effective  Homogeneous  Medium 


Figure  19:  Conceptual  schematic  of  the  homogenization  approach  of  the  Self-Consistent  micromechanical  model  for  a  coated 
inclusion. 


In  general,  the  SC  model  can  be  well  understood  as  an  intelligently  weighted  rule  of  mixtures 
approximation  as  shown  in  Eq.  (12), 

Cef[,  =  CMkl  +  f1 

ijkl  ijkl  J 

In  this  expression  the  stiffness  tensors,  Cx ,  of  each  material  are  denoted  with  an  M,  I,  or  C  superscript 

to  represent  the  matrix,  inclusion  or  coating,  respectively,  and  the  tensors  A7  and  Ac  represent 
quantities  known  as  strain  localization  tensors  for  the  inclusion  and  coating  phases,  respectively.  The 

strain  localization  tensors  relate  the  macroscopic  strain,  £9. ,  to  the  average  inclusion  (or  coating)  strain, 
£\7 ,  through  the  relation 


c1.  -  CM 

ijmn  ijmn 


A1 

mnkl 


+  /c 


Cc  -CM 

ijmn  ijmn 


mnkl 


(12) 


e1.  =  A1, 


ijkl  kl 


(13) 


For  the  current  problem,  the  volume  fraction  of  coating  material  is  set  equal  to  zero,  fc  =  0,  leaving 

only  the  inclusion  fraction  of  inclusion,  f1,  to  influence  the  overall  properties.  The  strain  localization 
tensors  represent  the  intelligent  weighting  parameters  and  have  significant  influence  on  the  effective 
properties  calculated.  They  take  into  account  not  only  the  material  anisotropy  of  each  constituent 
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phase,  but  also  the  inclusion  form.  For  the  special  case  of  no  coating,  the  inclusion  strain  localization 
tensor  reduces  to. 


A!jkl  ~  +  Tijmn  (Cm){C!nnkl  ~  Cmnkl)  ’ 


(14) 


where  I..kl  =  I4  is  the  fourth  order  identity  tensor  and  T1  ^Ce// j  is  the  volume  average  of  the  volume 
integral  of  the  what  is  known  as  the  modified  Green's  tensor,  G°jf  (r  —  r  ') , 


v‘  v 


Here  V1  represents  the  volume  of  the  inclusion  and  the  Green's  tensor  calculates  the  displacement  in 
the  /-direction  at  the  point  r  due  to  a  unit  force  in  the  ./-direction  at  resulting  from  an  inhomogeneity 

located  at  point  r'.  Ge/f  (r  —  r  is  defined  as  the  tensor  that  satisfies  the  differential  equation 


Full  evaluation  of  Eq.  (15)  is  a  non-trivial  task  and  details  of  the  Fourier  transform  technique  for  to  find 
an  approximate  solution  can  be  found  in  reference  [9].  Further,  inspection  of  Eq.  (14)  shows  that  the  SC 
model  is  implicit  and  numerical  methods  are  required  for  its  implementation.  However,  it  is  known  to 
provide  very  accurate  solutions  for  a  wide  range  of  problems  [7],  so  its  use  has  been  adopted  for  this 
project. 

Using  the  results  of  the  FEA  presented  in  the  previous  sections  to  find  the  effective  stiffness  of  the 
mesoscopic  inclusion.  Figure  20  shows  the  effective  properties  calculated  for  a  composite  containing  2% 
by  volume  of  identically  oriented  microstructured  NS  inclusions.  The  matrix  material  was  assumed  to  be 

isotropic  and  the  real  part  of  the  plane  wave  modulus  in  the  Xi-direction  (C^)  was  calculated  for 
various  properties  of  the  matrix,  C'n  /  M  ,  were  calculated  based  on  the  results  above  where 
maximum  negative  plane-strain  bulk  modulus  values  were  observed.  The  matrix  was  assumed  to  have  a 
small  loss  factor,  rf1  =  0.05  and  its  Poisson's  ratio  was  set  at  uM  =  0.30 .  The  results  present  the 
overall  properties  normalized  by  the  properties  of  the  matrix  for  a  wide  range  of  possible  inclusion  to 
matrix  stiffnesses.  The  objective  in  making  these  plots  is  to  determine  what  ratio  of  negative  inclusion 
stiffness  to  matrix  stiffness  provides  the  desired  overall  performance.  Inspection  of  the  resulting 
effective  anisotropic  properties  clearly  show  that  drastic  changes  in  overall  stiffness  and  loss  factor 

occur  for  stiffness  ratios  of  /  M  e  —0.25  —0.125  (approximately).  Note  that  the  scale  of 

the  normalized  loss  factor  is  logarithmic  in  order  to  clearly  display  the  large  range  of  increase  due  to  the 
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inclusion  of  negative  stiffness  particles  in  the  lossy  matrix.  The  largest  predicted  increase  occurs  for  the 
and  Cf!  —  ,  which  show  maximum  loss  factor  ratios  of  «  A0r]M  and  «  lOrf1 ,  respectively. 

It  s  worth  noting  that  these  plots  indicate  that  it  is  possible  to  simultaneously  increase  both  the  stiffness 
and  loss  of  a  composite  compared  to  the  matrix  material  by  using  just  a  very  small  amount  of 
constrained  NS  inclusions.  These  encouraging  results  suggest  that  viscoelastic  materials  containing  even 
very  low  volume  fractions  of  inclusions  structured  as  shown  in  Figure  10  will  display  drastic  increases  in 
energy  absorption  while  the  change  in  overall  stiffness  can  be  tuned  bases  on  the  ratio  of  the  inclusion 
to  matrix  stiffness. 


Figure  20:  Effective  stiffness  and  loss  properties  predicted  by  the  Self-Consistent  model  for  a  composite  containing  identically 
oriented  2%  by  volume  of  the  mesoscale  inclusions  modeled  in  Sections  4.1.5  and  4.1.6.  The  matrix  is  assumed  to  have  v  = 
0.30  and  q  =  0.05. 


4.1.8  Validation  of  the  Mesoscale  Effective  Stiffness  Approximation 

The  multi-scale  modeling  approach  developed  by  this  research  is  a  two-step  nonlinear  homogenization 
method.  It  employs  finite  element  models  to  simulate  the  force-displacement  behavior  of  the  structured 
inclusion  for  a  series  of  boundary  conditions.  The  FE  model  considers  the  nonlinear  elastic  response 
resulting  from  the  geometric  nonlinearity  on  the  inclusion  structure.  The  methodology  then  assumes 
that  the  nonlinear  stress  and  strain  behavior  resulting  from  the  inclusion  structure  can  be  well- 
represented  as  a  continuous  elastic  solid  inclusion  with  nonlinear  mesoscopic  effective  elastic 

properties,  C^so(Emeso),  where  C  is  the  effective  nonlinear  mesoscopic  stiffness  tensor  and  E  is  the 

Green's  strain  tensor  evaluated  at  the  mesoscale  .  This  is  the  micro-  to  meso  scale  transition  model.  The 
mesoscopic  stiffness  tensor  is  then  used  as  a  strain-dependent  input  to  the  meso-to-macro  scale 
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transition  model  which  is  based  on  standard  self-consistent  micromechanical  effective  medium  theory 
[7,8].  Using  this  approach,  one  is  thus  capable  of  direct  evaluation  of  the  influence  of  nonlinear 
microscale  structural  response  on  the  macroscopic  effective  behavior  of  a  composite  containing  a 
specified  volume  fraction  of  structured  inclusions.  To  validate  and  benchmark  this  unique  multiscale 
modeling  approach,  the  following  three  different  models  were  developed  for  determining  the  effective 
macroscopic  properties,  C^cro ,  of  a  composite  consisting  of  a  matrix  containing  with  a  known  volume 
fraction  of  structured  inclusions: 

1.  Micromechanical  effective  medium  models  of  a  matrix  containing  homogeneous  linear  elastic 
particulate  inclusions 

a.  The  three  micromechanical  models  used  were  the  Differential  Effective  Medium  (DEM), 
the  Self-Consistent  (SC),  and  the  Mori-Tanaka  (MT)  models. 

2.  FEA  of  representative  volume  element  (RVE)  consisting  of  a  matrix  containing  a  homogeneous 
cubic  inclusion  with  known  elastic  properties 

3.  FEA  of  a  RVE  consisting  of  a  matrix  containing  a  structured  inclusion 

The  figure  below  depicts  the  geometric  representation  of  these  3  model  sets. 


Figure  21:  Three  multiscale  modeling  validations  of  increasing  complexity,  (b)  Depicts  the  composition  and  geometry  used 
associated  with  analytical  modeling  called  for  in  bullet  1,  (c)  shows  the  geometry  and  properties  used  in  FEA  effective 
medium  models  associated  with  bullet  2,  and  (d)  depicts  the  geometry  and  properties  employed  doing  a  single-step 
multiscale  FEA  approximation  model  as  discussed  in  bullet  2. 
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Micro-  to  mesoscale  transition  homogenization  benchmarking  was  first  done  for  linear  inclusion 
behavior  and  was  validated  for  the  case  of  nonlinear  inclusion  behavior  by  extension.  In  order  to 
determine  C^so  (Emeso)  — >  C^soof  the  inclusion  for  input  into  model  sets  1  and  2,  a  linear  direct 

energy  method  was  carried  out  on  the  structured  inclusion  to  determine  the  small  strain  (linearized) 
stiffness  tensor  of  the  structured  inclusion.  The  effective  mesoscale  stiffness  results  for  the  suite  of 
applied  displacements  at  the  boundaries  of  the  inclusion  per  the  description  provided  in  Table  1  are 
given  in  Table  2. 


Table  2:  Strain  energy  and  stiffness  in  a  candidate  structured  inclusion  using  the  direct  FEA  micro-  to  mesoscale  transition 
model. 


Cn 

Cl2 

Cl3 

^33 

C44 

^66 

Strain  Energy  Function  (J) 

3.28E-10 

6.68E-10 

3.96E-06 

3.96E-06 

2.54E-07 

3.49E-08 

Stiffness  (Pa) 

1.87E+07 

3.44E+05 

4.28E+06 

2.26E+11 

3.62E+09 

4.98E+08 

The  meso-  to  macroscale  effective  medium  scale  transition  modeling  was  run  for  the  case  where  the 
matrix  material  was  considered  to  be  isotropic  constraining  with  Young's  modulus,  E  =  1  GPa  and 
Poisson's  ratio,  v  =  0.3.  Several  compositions  were  considered  with  volume  fractions,  cj),  of  structured 


inclusions  taking  values  of:  (j)  e  0.1  1,  2,  3, 


10 


% .  The  results  from  each  of  the  models  are 


shown  in  Figure  22-Figure  24. 


Effective  Axial  Stiffness  vs.  Inclusion  Volume  Fraction 


Figure  22.  Effective  axial  stiffnesses  versus  inclusion  volume  fraction.  Black  =  differential  effective  medium;  green  = 
generalized  self-consistent  model;  dark  blue  =  Mori-Tananka;  red  =  FEA  with  C7  =  C^so;  and  light  blue  =  FEA  with 

structured  inclusion  with  C7 . 
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Effective  Off-Diagonal  Stiffness  vs.  Inclusion  Volume  Fraction 


Figure  23.  Effective  off-diagonal  stiffnesses  versus  inclusion  volume  fraction.  Black  =  differential  effective  medium;  green  = 
generalized  self-consistent  model;  dark  blue  =  Mori-Tananka;  red  =  FEA  with  C7  =  ;  and  light  blue  =  FEA  with 

structured  inclusion  with  C7 . 

Effective  Shear  Stiffness  vs.  Inclusion  Volume  Fraction 


Figure  24.  Effective  shear  stiffnesses  versus  inclusion  volume  fraction.  Black  =  differential  effective  medium;  green  = 
generalized  self-consistent  model;  dark  blue  =  Mori-Tananka;  red  =  FEA  with  C7  =  C^so;  and  light  blue  =  FEA  with 

structured  inclusion  with  C7 . 

These  results  show  that  there  is  good  agreement  between  the  various  modeling  approaches  for 
determining  effective  macroscopic  stiffness  properties.  It  is  worth  noting  that  the  three  analytical  EMT 
models  produce  essentially  the  same  values  for  all  inclusion  volume  fractions  and  that  both  sets  of  FE 
model  results  closely  follow  these  trends.  At  the  extreme  volume  fraction  of  10%,  the  FEA-structured 
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inclusion  case  differs  from  the  EMT  results  by  no  more  than  20%.  For  the  purposes  of  this  research, 
where  volume  fractions  of  interest  are  limited  approximately  2%,  the  results  of  the  EMT  and  FEA 
approaches  are  much  more  similar.  We  can  therefore  be  certain  that  the  micro-  to  mesoscale  transition 
approach  employed  in  this  research  captures  the  physics  dominating  the  influence  of  microscale 
structure  on  the  effective  properties  at  the  macroscale.  Given  the  good  agreement  between  stiffness 
values  and  general  trending  among  all  models,  the  assumption  that  the  stress/strain  state  of  a 
structured  inclusion  can  be  well-represented  by  a  solid  inclusion  of  bulk  properties  is  thus  considered 
validated  and  of  utility  for  multiscale  design  purposes,  while  it  is  acknowledged  that  there  exists  levels  of 
discrepancy  that  are  higher  than  desired,  specifically  with  respect  to  the  effective  shear  moduli  of  the 
composite. 

4.1.9  Microstructural  Design:  Geometry  Parameterization 

In  light  of  the  ultimate  aim  of  designing  a  microscale  inclusion  tuned  for  desired  macroscale 
stiffness  and  loss  characteristics,  the  inclusion  geometry  was  parameterized  to  allow  for  a 
comprehensive  study  of  the  attainable  stiffness  values  given  the  material  and  thermal  loading 
constraints.  The  figure  below  illustrates  the  2D  inclusion  geometry  along  with  the  geometric  parameters 
of  interest.  The  external  boundaries  of  the  inclusion,  and  the  mesoscale,  is  defined  by  the  parameter  L . 
The  parameter  2 H  defines  the  height  of  T-shaped  interface.  The  parameter  T  defines  the  width  of  the 
connection  between  the  interface  and  the  buckled  element.  The  parameter  LI  is  defined  by  the  overall 
length  L\  =  L/2  —  2H  and  locates  the  midpoint  of  the  inclusion.  The  parameter  B  defines  the  height  of 
the  buckled  element  and  the  voided  region  below  it.  Finally,  L2  defines  the  buckled  element  length. 


Figure  25.  Cross-section  view  of  parameterized  inclusion  geometry  with  all  critical  parameters  labeled. 

Flaving  thus  established  the  various  geometric  parameters,  a  wide  range  of  inclusion  geometries  can  be 
constructed  by  means  of  implementing  a  select  number  of  parametric  ratios.  The  ratios  of  interest  are 


L2  2  B 

1112  LI'  ‘  L1-L2  ‘ 


and  R,~,t  =  —  . 

L2T  L2 


(17) 
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In  the  following  section,  these  ratios  are  varied  over  a  wide  range  of  possible  values  to  generate  a  highly 
populated  design  space  consisting  of  a  large  set  of  possible  geometries  that  may  produce  NS  behavior. 

4.1.10  Microstructural  Design:  Geometry  Dependent  Stiffness 

In  order  to  provide  a  full  inspection  of  the  design  space  for  this  particular  inclusion  shape,  a  range 
of  values  for  each  geometric  ratio  was  established  and  particular  geometries  were  constructed  using  a 
Halton  sequence.  A  Halton  sequence  is  a  deterministic  sequencing  approach  based  on  a  prime  number 
base.  The  Halton  sequence  offers  a  means  to  produce  quasi-random  sampling  of  a  wide  design  space, 
and  allows  one  to  easily  add  permutations  to  the  dataset,  provided  that  the  ranges  of  parameters  to  be 
inspected  remain  constant.  For  this  study,  the  number  of  geometric  trials  was  set  at  1000  within  the 
specified  ranges  of  values  that  the  parameters  can  take.  The  ranges  for  each  ratio  specified  for  design 
space  exploration  are  contained  in  Table  2.  These  ranges  were  chosen  to  provide  a  broad  sampling 
range  resulting  in  both  negative  and  positive  values  of  effective  mesoscopic  stiffness. 


Table  3:  Range  of  possible  values  taken  by  the  ratios  defined  in  Eq.  (17)  for  NS  inclusion  design  space  exploration. 


Ratio 

Min 

Max 

^LIL2 

0.6 

0.99 

rb 

0.1 

0.98 

R L2T 

0.05 

0.2 

Below  is  a  plot  of  each  parameter  value  for  any  given  trial.  The  data  points  illustrate  the  quasi-random 
nature  of  the  Halton  sequence  and  illustrate  the  degree  of  coverage  of  the  design  space  for  inspection. 
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Figure  26.  Parameter  value  versus  trial  number  for  the  ratios  used  to  generate  inclusion  geometries. 
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Using  Comsol  Livelink  for  MATLAB,  two  different  Finite  Element  Analysis  (FEA)  studies  were 
implemented.  Each  study  contained  an  automated  inclusion  geometry  construction  and  meshing 
routine.  The  first  study  analyzed  the  nonlinear  2D  plane  strain  temperature  dependent  elasticity 
constants  Cn,  C12,  and  C*66  over  a  range  of  imposed  displacements.  The  second  study  analyzed  the  3D 
elasticity  constants,  Cn,  C12,  C13,  C33,  C44,  and  ^66-  Below  are  nine  example  geometric  configurations  with 
associated  FEA  meshes  which  illustrate  the  wide  range  of  geometries  simulated  during  this  exhaustive 
design  space  inspection. 


Figure  27.  Example  geometries  generated  using  a  Halton  sequence  from  the  1000  trials  run  for  design  space  exploration. 


Below  is  a  plot  of  the  nonlinear  plane  strain  stiffness  values  (Cn,  Cl2/  and  C66)  versus  applied 
displacement  for  all  1000  geometric  simulations. 
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Figure  28.  Stiffness  versus  displacement  values  for  Cn,  C12,  and  C66  for  all  1000  trials. 


A  zoomed  section  of  the  displacement  dependent  mesoscopic  Cn  verifies  that  negative  values  for  the 
Cn  stiffness  were  obtained  for  a  percentage  of  the  designs  generated  during  the  design  space 
exploration. 
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Figure  29.  Zoomed  view  of  Cn  stiffness  versus  displacement  values  for  each  of  the  1000  trials  run  during  this  research. 

This  data,  along  with  values  of  the  3D  elastic  constants,  which  altogether  amounted  to  over  500,000 
data  points,  was  used  to  build  a  surrogate  model  of  the  stiffness  data  as  a  function  of  the  structured 
inclusion  parameters.  A  surrogate  model  is  an  analytical  function  that  is  fitted  to  arbitrary  data  points,  in 
this  case  points  generated  using  computationally  expensive  FEA  models,  in  order  to  save  time  and 
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computational  expense.  There  are  many  approaches  to  generating  surrogate  models  and  the  approach 
used  in  this  project  is  known  as  "Kriging"  [9a]. 

4.1.11  Microstructural  Design:  Kriging  Surrogate  Modeling 

Kriging  surrogate  modeling  is  a  data  fitting  technique  that  was  initially  developed  in  the  field  of 
topology  to  interpolate  the  value  of  a  random  field  at  an  unobserved  location  using  observations  of  its 
value  at  nearby  locations.  It  can  be  simply  described  as  a  known  global  function  with  additional 
departures  from  the  base  model,  as  described  in  Eq.  (18) 

y(x)=T;Pifi(x)+z(x)-  (18) 

7=1 

In  this  expression,  /?  are  unknown  coefficients  of  the  base  fitting  function,  f.  (x)  are  pre-defined 
functions,  and  Z (x)  provides  the  departures  from  the  base  function  in  order  to  interpolate  the  training 
points  [9b].  Z(x)is  the  realization  of  a  stochastic  process  with  a  mean  of  zero,  variance  of  <r2,  and 
nonzero  covariance  of  the  form: 


cov  Z(xi),z{xj)  =(T2R(xi,xj) , 


(19) 


where  R  is  a  user-specified  correlation  function.  Similar  to  the  approach  employed  by  Backlund  et  al. 
[9b],  this  study  employs  a  constant  term  for  /(x)  and  a  Gaussian  curve  in  the  form  shown  in  (20)  is 
employed  as  the  correlation  function: 


/  \ 

r  _  /  x2i 

( x. ,  Xj  j  —  exp 

- 1 

1 

1 

_ 1 

(20) 


In  Eq.  (20)  Ot  are  unknown  correlation  parameters  that  are  automatically  determined  during  the  model 
fitting  process. 

While  Kriging  has  been  shown  to  be  slow  in  terms  of  build  and  prediction  time  compared  with 
other  surrogate  modeling  techniques  [9c],  Kim  et  al.  demonstrate  the  ability  of  Kriging  to  successfully 
model  nonlinear  and  multimodal  problems  [9d].  For  this  reason,  a  Kriging  approach  was  selected  for 
generating  a  surrogate  model  of  the  nonlinear  effective  mesoscopic  inclusion  stiffness  as  a  function  of 
the  geometric  parameters  of  interest  to  this  project. 

In  general  terms,  surrogate  models  were  developed  for  this  project  using  a  Kriging  approach  in 
order  to  determine  the  effective  nonlinear  mesoscopic  stiffness  as  a  function  of  inclusion  geometry, 
constituent  material  properties,  and  thermo-mechanical  loading.  This  can  be  stated  in  a  pseudo- 
mathematical  representation  as  illustrated  below: 
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CmLo  =  /  ( geometry,  materials,  AT,  displacement ) . 


The  advantage  of  this  functional  approach  is  that  the  effective  mesoscopic  stiffness  can  be  calculated 
using  an  analytical  function  rather  than  a  computationally  expensive  FE  model.  Clearly,  this  is  a  powerful 
tool.  It  permits  the  estimation  of  the  mesoscale  inclusion  stiffness  for  any  arbitrary  geometry  (within  the 
provided  bounds)  in  a  matter  of  seconds  as  opposed  to  the  hours  required  to  generate  results  using 
computationally-intensive  FEA  methods.  It  is  important  to  note,  however  that  one  must  generate  a  large 
dataset  using  FEA  results  in  order  to  produce  the  surrogate  model.  Further,  the  models  are  established 
in  order  to  generate  the  training  data  points  for  subsequent  surrogate  models  should  inclusion  materials 
or  temperatures  change.  Thus,  a  very  general  approach  to  obtaining  accurate  estimates  of  inclusion 
stiffness  based  on  arbitrary  geometries  has  been  established  and  leveraged  in  this  work  to  generate 
significant  computational  capabilities  for  a  highly  nonlinear  constitutive  microscale  geometry  that  can 
be  exploited  to  produce  negative  stiffness  behavior.  This  will  greatly  enhance  the  future  design  of 
multiscale  materials  that  employ  microscale  structure  to  generate  negative  stiffness  behavior. 

4.2  Negative  Poisson  Ratio  Materials 

In  recent  years  there  has  been  an  increased  general  interest  in  a  group  of  materials,  the  properties 
of  which  have  been  described  as  auxetic  [10].  An  auxetic  material  is  one  in  which  the  Poisson's  ratio,  v,  is 
negative,  i.e.  an  applied  tensile  strain  in  the  longitudinal  direction  results  in  a  positive  strain  in  the 
transverse  direction.  It  has  been  found  that  a  number  of  naturally  occurring  materials  are  auxetic  in 
behavior,  for  example,  a  range  of  crystalline  structures  including  iron  pyrites  [11],  [12]  and  some  forms 
of  skin  [13]  and  cancellous  bone  structures  [14],  though  it  should  be  noted  that  Poisson's  ratios  for  bone 
are  usually  reported  as  positive  [15].  However,  it  is  also  possible  to  engineer  auxetic  behavior  and  this 
has  been  shown  to  be  possible  for  a  number  of  synthetic  materials.  These  have  included  polyurethane 
and  metallic  foams  [16-21],  composite  laminates  [22]  and  a  range  of  microporous  polymers  namely 
polytetraflouroethylene  (PTFE)  [23],  ultrahigh  molecular  weight  polyethylene  (UHMWPE)  [24],  nylon 
[25],  polypropylene  (PP)  [26],  and  polyester  [27]. 
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Figure  30:  Auxetic  materials  and  structures  across  the  length  scales. 


Considering  polymers,  it  was  found  that,  for  larger  cross  section  cylinders,  the  auxetic  effect  is  due 
to  their  complex  microstructure  [23-26].  This  is  essentially  a  network  structure,  which  has  been 
described  as  a  system  of  nodules  interconnected  by  fibrils  and  the  mechanism  for  this  effect  is  well 
defined  [28],  [29].  In  essence,  the  re-entrant  structure  of  auxetic  foams  had  been  reproduced  by  the 
nodule-fibril  structure  [26],  the  schematic  structure  for  which  is  shown  in  Figure  31  and  a  similar 
deformation  mechanism  was  noted.  In  this  case,  it  is  the  fibrils  which  are  hinging,  though,  rather  than 
the  cell  ribs  as  is  seen  in  the  foams.  This  hinging  of  the  fibrils  in  turn  causes  the  nodes  to  translate, 
producing  the  auxetic  effect.  It  was  also  noted  that  there  was  a  variation  in  the  number  and  size  of  the 
fibrils  formed  during  the  production  of  the  cylinders  for  the  various  polymers,  in  particular,  the  degree 
of  fibrillation  for  PP  was  much  lower  than  that  for  PTFE  or  UFIMWPE  [26]. 
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Figure  31:  Two  dimensional  schematic  diagram  of  nodule-fibril  model  for  auxetic  behavior  in  polymeric  cylinders 


4.2.1  Processing  routes  to  produce  auxetic  extruded  fibers  and  films 

In  2001,  a  novel  processing  route  was  developed  to  produce  auxetic  polypropylene  fibers  based  on 
the  powder  processing  route  used  to  produce  auxetic  PP  cylinders  [30], [31].  The  same  powder  and 
processing  temperature  were  employed,  i.e.  Coathylene  PB0580  and  at  a  temperature  of  159°C.  The 
powder  processing  technique  was  adapted  to  take  place  within  an  Emerson  and  Renwick  Ltd  Labline 
melt  extruder.  The  extruder  has  5  temperature  zones,  each  with  individual  thermostatic  control,  a  3:1 
compression  ratio,  a  25.4mm  (1  inch)  diameter  screw  and  a  length/diameter  ratio  of  24:1.  A  schematic 
of  the  extruder  is  shown  in  Figure  32.  Powder  is  fed  into  the  hopper  and  is  transported  through  the 
extruder  by  the  screw.  In  conventional  melt  extrusion,  the  polymer  is  fully  molten.  In  this  case,  however, 
processing  occurs  at  a  much  lower  temperature  with  the  aim  being  to  partially  melt  the  powder  only.  To 
produce  fibers,  the  extruder  is  fitted  with  a  spinneret  plate.  The  first  fibers  were  produced  using  a  40 
hole  spinneret  plate  with  hole  diameters  of  approximately  0.55mm  [30], [31].  The  fibers  produced  were 
characterized  by  measuring  their  Poisson's  ratio  using  Biax  200C  micro-tensile  testing  equipment,  which 
was  fitted  with  a  10  Newton  load  cell  and  specialized  Messphysik  videoextensometry  software,  shown 
schematically  in  Figure  33. 
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Figure  33:  Schematic  of  the  videoextensometry  measurement  system  used  to  characterize  auxetic  fibers. 


The  Messphysik  software  simultaneously  measures  length  and  width  data  from  changes  in  contrast 
between  attached  markers  (along  the  length)  and  the  edges  (across  the  width)  of  the  fibers.  This 
allowed  the  optimum  conditions  to  be  defined  as  a  processing  temperature  of  159°C,  screw  speed  of 
1.047rads-1  (equivalent  to  10  RPM  and  resulting  in  a  throughput  of  6g/mm)  and  a  take-off  speed  of 
2m/min.  No  drawing  other  than  that  caused  by  gravity  as  the  fiber  exited  the  spinneret  was  employed. 
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Using  these  conditions,  40%  of  the  fibers  tested  displayed  auxetic  character.  These  fibers  were  found  to 
be  inhomogeneous  along  their  length  with,  on  average,  20  -  30%  of  their  lengths  consisting  of  auxetic 
material  i.e.  their  auxeticity  is  20  -  30%.  The  Poisson's  ratio  of  the  auxetic  region  was  found  to  be,  on 
average,  v  =  -1.62  and  the  overall  Poisson's  ratio  of  the  entire  fiber  length  was  also  negative,  being  v  =  - 
0.7  [32].  The  auxetic  fibers  produced  are  undrawn,  and,  as  such,  as  well  as  being  inhomogeneous,  they 
have  low  mechanical  properties,  with  a  modulus  of  1.3  GPa,  strain  to  yield  of  7.5%  and  tensile  strength 
of  25MPa 

For  the  fibers,  auxetic  behavior  is  based  on  a  closely  packed  particulate  or  granular  microstructure 
with  a  low  void  content,  which  is  a  close  approximation  to  the  actual  microstructure,  see  Figure  34. 


Figure  34:  Microstructure  of  an  auxetic  polypropylene  fiber. 


The  nodule-fibril  model  used  to  study  auxetic  behavior  in  cylinders  is  clearly  not  appropriate  here.  In  this 
case,  mechanisms  for  auxetic  behavior  are  based  on  a  rough  particle  model  described  by  interlocking 
rigid  hexagons  as  shown  in  Figure  35  [33],  [34],  The  interlocked  structure  is  formed  by  rectangular  radial 
'keys'  and  geometrically  matched  'keyways'  and  is  derived  from  the  geometry  employed  in  the  Tokai 
nuclear  reactor  designed  in  the  1950's  [35],  [36].  These  structures  consist  of  free-standing  columns  of 
graphite  bricks  laterally  connected  by  loose  side  and/or  corner  keys  in  keyways.  The  key-keyway 
combinations  produce,  effectively,  a  rough  macroparticle  assembly  in  which  interactions  are 
characterized  by  having  tangential  stiffness  exceeding  normal  contact  stiffness.  The  keyed-brick 
structures  expand  in  all  radial  directions  when  subjected  to  a  tensile  load  in  the  horizontal  plane  by 
translation  of  the  bricks  through  sliding  of  the  bricks  and  radial  keys  along  the  keyways  [34].  To  study 
the  auxetic  fibers,  the  interlocking  hexagon  geometry  has  been  selected  for  simplicity. 


Final  Report 


New  Solutions  for  Energy  Absorbing  Materials 


Page  42  of  68 


Figure  35:  Schematic  of  the  interlocking  hexagon  model.  In  (a),  the  material  is  fully  densified;  in  (b),  the  material  is  partially 
expanded. 


Looking  at  Figure  35,  the  hexagons  have  edge  lengths  and/2 .  The  edges  of  are  aligned  parallel 
to  the  x-axis  and  the  edges  of  l  are  at  an  angle  a  to  the  x-axis.  Thus,  the  unit  cell  lengths  are  given  by 


and 


X  =  ik+l,_  cos  a  +  aj , 
Y  =  2  (l2  sin  a  +  a  cot  a ) , 


(21) 

(22) 


where  a  is  the  parameter  defined  in  Figure  35  as  the  gap  between  adjacent  hexagon  particles. 

For  deformation  of  the  structure  by  the  mechanism  of  particle  translation,  a  is  the  variable 
parameter.  So,  the  changes  in  unit  cell  lengths  are  thus 

dX  n  dY  n  x 

- =  2  and - =  2  cot  a ,  (23) 

da  da 

and  the  Poisson's  ratio  in  the  x-direction  is  given  by 
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-  de«  dY_[Y_  [dY]{dX]  ’  [x 

Vxy~  de  dX  /  X~  da  da  Y 

Combining  Eqs.  (21)-(24)  one  finds  ^  to  be 

xy 

cos  a  U  +  l2  cos  a  +  a) 

vx  = - — - - - , 

xy  l2  sin2  a  +  a  cos  a 

and  similarly  u  is 

yx 

l  sin2  a  +  a  cos  a 

V  x  = - — - r  . 

cos  aU  +l2  cos  a  +  a) 


(24) 


(25) 


(26) 


The  expressions  for  Young's  moduli  of  this  network  were  derived  by  using  the  conservation  of 
energy  approach  [34].  In  the  x-direction,  £x  is  given  by 

\  / 

2  cos2  a,  +  1  ^  +  l2  cos  ol  +  a 

sin  a  l  sin2  a  +  a  cos  a 

/  \  z 

where  kh  is  the  stiffness  of  a  spring  assumed  to  connect  each  male-female  interlock  combination. 
Similarly,  the  Young's  modulus  in  the  y-direction  is 

,  \(i  ■  2  t 

„  ,  2  cos  a  +  1  L  sin  a  +  a  cos  a 

E  =k  - —  ^ -  .  (28) 

v  1  sin  a  cos2  a  l  +  l  cos  a  +  a 

V  J  \  1  2  y 

Predictions  show  that  interlocking  hexagons  having  0<a<90°  lead  to  auxetic  behavior  for  the  assembly. 
The  models  show  good  agreement  with  experimental  values  at  low  strains  [33],  [34]. 

A  further  experimental  variation  on  the  extrusion  process  was  introduced  by  changing  the 
spinneret  plate  employed  to  produce  fibers  to  a  slit  orifice  of  dimensions  63.5  x  14.2  x  0.38mm.  Using 
this,  films  were  produced  from  polypropylene,  initially  based  on  the  process  window  that  was  used 
successfully  to  produce  PP  fibers  [27],  [33].  The  PP  films  produced  at  159°C  with  screw  speed  1.05  rads'1 

and  0.0225  ms'1  take  off  speed  were  found  to  have  =  -1.12  and  u  =  -0.77  at  low  strain,  typically 

up  to  1%.  The  Young's  moduli  in  this  case  were  £x  =  0.34  GPa  and  £y  =  0.20  GPa.  Beyond  1%,  the  film 
undergoes  an  auxetic-to-nonauxetic  transition  [27],  A  detailed  processing  parameter  study  was 
undertaken  to  assess  the  effect  of  the  change  in  cross-section  of  the  die-head  from  circular  for  the  fibers 
to  rectangular  for  the  films  [37],  There  were  several  interesting  findings.  For  the  first  time,  a  100% 
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auxetic  extruded  product  was  produced  for  several  processing  parameter  sets.  For  example,  100% 
auxeticity  was  measured  for  films  extruded  at  159°C,  screw  speed  1.05  rads 1  and  take  up  speed  0.0225 
ms'1  and  for  those  extruded  at  159°C,  screw  speed  1.05  rads 1  and  take  up  speed  0.15  ms'1. 

4.2.2  Information  on  NPR  materials  selected  and  why  they  are  of  interest 

The  inclusions  selected  for  this  project  are  a-cristobalite  and  samples  produced  from  auxetic 
extruded  products  comprising  polypropylene  (PP)  fibers  or  films.  The  a-cristobalite  was  selected  as  it 
has  a  high  modulus,  E  =  65.4  GPa,  combined  with  a  small  negative  Poisson's  ratio  of  v  =  -0.16. 
Conversely,  the  modulus  of  the  auxetic  fibers  and  films  is  much  lower  at  E  =  1.3  GPa  for  the  fibers  and  E 
=  0.34  GPa  for  the  films,  but  the  Poisson's  ratio  in  both  cases  is  higher.  On  average  for  both  systems  the 
Poisson's  ratio,  which  is  strain  dependent,  is  v  =  -0.6.  The  fibers  can  have  a  Poisson's  ratio  as  low  as  v  =  - 
1.6  with  auxeticity  of  20  -  30%.  The  films,  while  of  a  lower  modulus,  can  be  produced  to  have  v  =  -0.95 
and  be  100%  auxetic.  Thus,  the  inclusions  selected  for  this  project  are  available  or  can  be  readily 
manufactured,  and  comprise  a  high  modulus,  low  negative  Poisson's  ratio  inclusion  or  a  low  modulus, 
high  negative  Poisson's  ratio  inclusion. 

4.2.3  Fabrication  of  a-cristobalite  samples 

Samples  have  been  manufactured  to  the  specification  of  0,  5,  10,  15,  20  and  25wt%  a-cristobalite, 
which  is  a  naturally  occurring  high  modulus  auxetic  mineral.  The  resin  used  was  a  room  temperature 
curing  epoxy  system  (Epon  E828/D400)  together  with  a  viscosity  modifier  (silicone  polycarbonate 
urethane,  CarbiSil®)  in  trace  amounts.  The  inclusions  were  added  to  the  resin,  vigorous  mixing  was 
undertaken  for  1  hour  and  then  the  homogeneous  slurry  was  cast  in  aluminum  moulds  and  cured  at 
room  temperature. 

4.2.4  Microscopy 

Optical  microscopy  was  not  successful  in  revealing  the  microstructure  of  the  samples  due  to  their 
opacity,  as  can  be  seen  in  Figure  36  and  Figure  37  which  optical  micrographs  of  the  0wt%  and  25wt% 
samples,  respectively. 
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Oft  a  -  C  rysta  barite  tap  palsihed  surface 


Q 94  a  -  Crys ta  b aFrte  tap  scratched  surface 


Figure  36:  Optical  micrograph  of  a  0%  a-cristobalite  sample. 


2  5  ^  a  -  Csrys  tail  a  t'rte  tap  pafished  surface 


2594a-Qrystat]afte  side  surface 


259i  a  -  Cjysta  hafte  ta  p  scratched  s  urface 


259i  a  -  Ctysta  baCte  scratched  side  surface 


Figure  37:  Optical  micrograph  of  a  25%  a-cristobalite  sample. 
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It  was  therefore  decided  that  scanning  electron  microscopy  be  carried  out  in  place  of  optical  microscopy 
and  that  this  should  be  in  conjunction  with  elemental  mapping  to  begin  a  study  of  the  distribution  and 
interfacial  bonding  of  the  a-cristobalite  particles.  The  findings  of  the  elemental  mapping  are  given  in 
Table  4  and 
Table  5. 


Table  4:  Weight  %  of  carbon,  oxygen  and  silicon  in  the  a-cristobalite  inclusion  samples 


%  a-cristobalite 

Weight  %  of 

C 

O 

Si 

0 

75.18 

23.69 

- 

5 

74.21 

24.35 

1.05 

10 

71.85 

26.55 

1.61 

15 

70.405 

26.605 

3.09 

20 

70.69 

25.54 

3.77 

25 

68.21 

26.62 

5.03 

Table  5:  Atomic  %  of  carbon,  oxygen  and  silicon  in  the  a-cristobalite  inclusion  samples. 


%  a-cristobalite 

Atomic  %  of 

C 

O 

Si 

0 

80.74 

19.10 

- 

5 

79.60 

19.93 

0.48 

10 

77.71 

21.55 

0.74 

15 

76.84 

21.72 

1.45 

20 

77.28 

20.97 

1.76 

25 

75.46 

22.11 

2.38 

Table  6  and  Table  7  show  that  while  the  distribution  is  not  homogeneous  throughout,  the  particles  do 
not  float  to  the  top  of  the  samples  or  sink  to  the  bottom,  creating  intensely  rich  a-cristobalite  inclusion 
regions.  The  data  presented  in  these  tables  indicate  that  SEM  micrographs  allow  for  a  clearer  picture  of 
the  internal  structure  of  the  samples. 


Table  6:  Weight  %  of  carbon,  oxygen  and  silicon  in  the  top,  bottom  and  middle  of  the  25%  a-cristobalite  inclusion  samples 


Surface  examined 

Weight  %  of 

C 

O 

Si 

Top 

71.72 

26.93 

1.27 

Middle 

68.21 

26.62 

5.03 

Bottom 

71.31 

26.89 

1.80 
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Table  7:  Atomic  %  of  carbon,  oxygen  and  silicon  in  the  top,  bottom  and  middle  of  the  25%  a-cristobalite  inclusion  samples 


Surface  examined 

Atomic  %  of 

C 

0 

Si 

Top 

77.52 

21.89 

0.59 

Middle 

75.46 

22.11 

2.38 

Bottom 

77.29 

21.88 

0.84 

To  illustrate  the  structures  seen  during  SEM  imaging  of  the  samples  constructed  for  this  project,  two 
micrographs  are  presented  below.  The  first,  Figure  38,  shows  the  side  view  of  a  resin  only  sample. 


Figure  38:  SEM  image  of  a  pure  resin  sample  for  reference  with  composite  consisting  of  resin  and  a-cristobalite  inclusions 
shown  in  Figure  39. 

The  second  image  shows  the  side  view  of  a  25wt%  composite.  The  red  circles  indicate  the  locations  of  a- 
cristobalite  inclusions,  confirmed  by  the  elemental  mapping.  The  resulting  changes  in  material 
properties  for  various  samples  are  discussed  in  Section  4.2.5. 
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Figure  39:  SEM  image  of  a  composite  consisting  of  resin  and  25wt%  a-cristobalite  inclusions.  Red  circles  indicate  inclusions. 


4.2.5  Quasistatic  Mechanical  Testing  of  a-cristobalite  Inclusion  Composites 

An  average  of  6  samples  were  tested  for  each  of  the  composites  manufactured  with  0,  5,  10, 15,  20 
and  25wt%  a-cristobalite,  each  specimen  being  150mm  long  by  15mm  wide  by  5mm  thick.  Tensile  tests 
were  carried  out  to  measure  the  mechanical  properties  of  the  composites  using  an  Instron  3369,  in 
conjunction  with  videoextensometry  using  a  Messphysik  ME  46  videoextensometer.  A  crosshead  speed 
of  0.3mm/min  was  used  with  a  gauge  length  of  100mm  and  a  load  cell  of  100N.  The  Poisson's  ratio 
measurement  was  obtained  using  a  rectangular  box  of  30  by  10mm  marked  on  the  samples  for  the 
videoextensometer  to  track.  The  strain-strain  plots  for  the  control  sample  (i.e.  with  no  a-cristobalite) 
and  for  the  25%  a-cristobalite  sample  are  shown  below  in  Figure  40.  In  both  cases,  the  Poisson's  ratio  is 
positive  but  it  can  be  seen  that  the  addition  of  the  25wt%  a-cristobalite  reduces  the  Poisson's  ratio 
value.  In  the  cases  shown,  this  reduction  is  from  v  =  +0.33  to  v  =  +0.20  for  the  a-cristobalite  inclusion 
sample.  A  complete  set  of  the  results  from  these  tests  is  given  in  Table  8. 

The  observed  reduction  in  Poisson's  ratio  is  in  agreement  with  the  trend  shown  by  theoretical 
predictions  carried  out  using  the  self-consistent  (SC)  analytical  model  and  FEA,  both  of  which  are  shown 
in  Figure  41  together  with  the  experimental  data.  The  system  used  here  has  £f/£m  =  24.  The  effect  of 
using  £f/£m  =  40  is  also  shown  in  the  same  graph  given  that  it  is  of  interest  since  it  should  cause  the 
Poisson's  ratio  to  be  further  reduced.  Agreement  is  very  good  at  low  volume  fraction  of  a-cristobalite. 

Similar  plots  were  obtained  for  the  Young's  modulus,  and  these  can  be  seen  in  Figure  42.  Here,  the 
experimental  data  are  compared  to  FE  analysis  and  analytical  values  calculated  using  the  SC  and  Hashin- 
Shtrikman  (HS)  upper  and  lower  bounds.  For  ease  of  comparison,  the  modulus  value  plotted  is  that  of 
the  composite  relative  to  the  matrix.  The  addition  of  a-cristobalite  leads  to  a  rise  in  the  Young's  modulus 
of  the  composite,  with  the  HS  lower  bound  predictions  agreeing  very  well  with  the  experimental  values. 
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True  Strain  in  X-direction 

(a) 


(b) 


Figure  40:  True  strain-strain  plot  for  the  (a)  control  and  (b)  25wt%  composites. 


Table  8:  Experimental  Poisson's  ratio  values  for  the  a-cristobalite  composites. 


Sample 

Specimen  number 

Average 

Poisson's 

ratio 

1 

2 

3 

4 

5 

6 

Resin 

0.44 

0.33 

0.24 

0.30 

0.20 

0.32 

0.32±0.09 

Resin  + 

5% 

0.33 

0.16 

0.18 

0.18 

0.33 

0.26 

0.24±0.08 

Resin  + 

10% 

0.22 

0.18 

0.19 

0.24 

0.28 

0.28 

0.23±0.04 

Resin  + 

15% 

0.20 

0.17 

0.17 

0.23 

0.25 

0.30 

0.22±0.05 

Resin  + 

20% 

0.11 

0.24 

0.26 

0.18 

0.21 

0.21 

0.20±0.05 

Resin  + 

25% 

0.10 

0.18 

0.17 

0.19 

0.28 

0.27 

0.20±0.07 
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Figure  41:  Plot  of  Poisson's  ratio  against  volume  fraction  of  a-cristobalite  inclusions  for  £f/£m  =  24  and  40. 
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Figure  42:  Plot  of  relative  Young's  modulus  against  volume  fraction  of  a-cristobalite  inclusions 


A  number  of  reasons  for  the  disagreement  between  the  modeling  and  experimental  results  have 
been  considered.  The  first  of  these  is  that  the  models  are  idealized.  The  second  is  that  we  think  that 
despite  the  use  of  the  viscosity  modifier,  the  particle  distribution  of  the  cristobalite  is  not  perfect, 
leading  to  some  particle  to  particle  contact,  especially  as  the  %  volume  of  the  cristobalite  increases.  The 
models  are  based  on  complete  distribution  with  no  particle  to  particle  contact,  which  would  tend  to 
reduce  the  properties.  The  interface  on  inspection  in  the  SEM  appears  to  be  a  further  possible  source  of 
difference  between  the  model  and  experimental  values.  The  project  was  stopped  prior  to  detailed 


Final  Report 


New  Solutions  for  Energy  Absorbing  Materials 


Page  51  of  68 


examination  of  these  hypothetical  problems,  but  results  clearly  suggest  that  the  samples  are  imperfect, 
which  is  not  unexpected.  An  idealized  bond  between  particle  and  matrix  is  an  assumption  of  the  models 
used  for  comparison,  so  the  interface  should  be  a  focus  of  future  work  on  these  types  of  composite 
materials  in  order  to  gain  maximum  properties.  A  final  reason  for  the  observed  difference  could  be  that 
the  samples  have  reduced  volume  fractions  than  those  reported.  The  samples  were  made  externally  to 
our  recipe  but  the  initial  elemental  analysis  appears  to  indicate  that  there  may  be  less  Si  present  than 
was  desired.  Eliminating  the  sources  of  all  of  these  minor  discrepancies  should  be  the  focus  of  future 
work  in  the  area  of  auxetic  composite  materials. 

4.2.6  Energy  absorptive  capacity  of  the  composites 

There  are  many  references  which  support  the  increased  energy  absorptive  capacity  of  auxetic 
materials  in,  for  example,  foams  [38], [39]  and  carbon  fiber  laminates  [22].  More  specifically  for  this 
project,  investigations  into  the  properties  of  auxetic  UHMWPE  cylinders  were  carried  out  and  they 
revealed  benefits  in  deploying  auxetic  polymeric  materials.  An  interesting  study  was  performed  to 
measure  the  dynamic  modulus  of  the  auxetic  UHMWPE  using  ultrasound  [40]  which  resulted  in 
measurement  of  the  attenuation  coefficient,  a,  for  auxetic,  conventional  compression  molded  plaques 
and  compacted  and  single  sintered  UHMWPE.  The  microporous  polymers,  whether  auxetic  or  sintered 
(which  resulted  in  a  positive  Poisson's  ratio),  showed  very  large  enhancements  in  the  attenuation 
coefficient,  a.  The  authors  report  purely  on  experimental  values,  however,  that  the  auxetic  form  was  so 
good  at  absorbing  the  ultrasonic  signal  that  in  the  highly  fibrillar  form,  it  was  not  possible  to  obtain  a 
single  measurement  of  a  as  the  signal  was  completely  absorbed  after  1  single  back  reflection. 

The  second  of  the  directly  related  properties  which  has  been  studied  in  detail  was  the  indentation 
resistance,  H.  This  is  related  to  the  Poisson's  ratio  as 


H  oc(l-u2)X ,  (29) 

with  x  dependent  on  the  method  of  analysis  considered.  For  example,  x  =  2/3  for  Hertzian  indentation. 
Ball  indentation  tests  were  carried  out  on  auxetic  UHMWPE  using  conventional  compression  molded 
plaques  and  compacted,  single  sintered  UHMWPE  as  comparison  material.  It  was  found  that  indentation 
resistance  was  enhanced  by  up  to  a  factor  of  3  at  low  loads,  where  the  strain  dependent  negative 
Poisson's  ratio  has  its  highest  values.  The  mechanisms  at  play  in  the  indentation  response  were  then 
examined  [41].  Alderson,  Fitzgerald  and  Evans  propose  that  the  fibrils  pull  the  nodules  under  the 
indenter,  leading  to  local  and  elastic  densification  under  the  indenter.  This  is  underlined  by  the 
observation  that  the  indents  in  the  auxetic  UHMWPE  recover  whereas  the  other  conventional  forms 
show  permanent  plastic  deformation.  The  recovery  is  suggested  to  be  due  to  the  fibrils  unwrapping  on 
removal  of  the  indenter,  allowing  the  nodules  to  return  to,  or  close  to,  their  original  position.  Though 
these  previous  studies  showed  increases  in  energy  absorption  when  auxetic  materials  were  tested,  it  is 
unclear  whether  particulate  composite  materials  consisting  of  a  conventional  matrix  material  containing 
auxetic  inclusions  will  also  display  increases  in  energy  absorption  capacity.  The  materials  that  introduced 
in  Section  4.2.3  and  tested  quasistatically  in  Section  4.2.5  were  further  characterized  using  dynamic 
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mechanical  analysis,  modal  damping,  and  ultrasonic  methods.  The  results  of  those  measurements  are 
given  in  the  following  section. 

4.2.7  Characterization  of  Particulate  Composites  Containing  a-critobalite 

The  auxetic  inclusion  composites  fabricated  by  the  University  of  Bolton  were  characterized  across  a 
broad  range  of  frequencies  to  capture  their  dynamic  stiffness  and  loss  performance.  This  was 
accomplished  though  Dynamic  Mechanical  Analysis  (DMA),  model  damping,  and  ultrasonic 
characterization  techniques.  These  tests  were  run  by  researchers  at  all  three  of  the  associated 
universities.  The  University  of  Bolton  characterized  the  samples  using  DMA,  the  University  of  Bristol 
measured  the  dynamic  response  using  forced  vibration  techniques,  and  the  University  of  Texas  at  Austin 
measured  the  response  using  contact  ultrasonic  techniques.  The  results  of  these  measurements  are 
provided  in  the  following  subsections. 

3.2.7.1  Dynamic  mechanical  analysis 

Mechanical  properties  of  materials  can  be  expressed  via  dynamics  modulus  represented  with  the  a 
storage  and  loss  modulus  of  the  material  [42].  One  conventional  means  for  measuring  these  quantities  is 
through  Dynamic  mechanical  analysis  technique.  Dynamic  Mechanical  Analysis  is  a  well-established, 
robust  technique,  which  is  frequently  used  to  characterize  the  viscoelastic  behavior  of  the  polymers  and 
polymeric  composite  materials  [43].  When  a  material  is  tested  using  the  DMA  technique,  an  oscillatory 
(sinusoidal)  strain  (or  stress)  is  applied  to  the  material  and  the  resulting  stress  (or  strain)  developed  in 
the  material  along  with  its  relative  phase  to  the  imposed  loading  is  measured  to  yield  a  complex 
modulus.  The  complex  modulus  is  thus  the  frequency-dependent  material  resistance  to  the 
deformation.  It  has  two  components,  the  storage  modulus,  which  is  a  parameter  describing  the  elastic 
energy  stored  during  a  single  cycle  of  loading,  and  the  loss  modulus,  which  is  a  parameter  of  the  energy 
dissipated  through  conversion  to  heat  in  a  single  loading  cycle.  The  Dynamic  Mechanical  Analyzer 
produces  very  important  information  on  all  changes  in  the  state  of  molecular  motion  as  temperature  or 
frequency  is  scanned  and  can  be  used  to  characterize  material  over  a  very  broad  range  of  temperatures 
and  frequencies. 

DMA  test  was  performed  at  the  University  of  Bolton  using  dual  cantilever  clamp  on  Isothermal 
frequency  sweep  mode  for  the  series  of  a-cristobalite  samples.  Because  the  stiffness  of  the  samples  is 
moderate,  dual  cantilever  clamp  was  used  for  the  experiment.  The  sample  sizes,  were  60mm  x  15mmx 
5mm,  were  determined  based  on  the  defined  standards  to  insure  accurate  measurements.  The 
isothermal  frequency  sweep  test  was  performed  at  35°C  under  15pm  amplitude  which  limited 
deformation  to  the  linear  regime.  The  storage  modulus  and  tan  6  (ratio  of  the  loss  to  storage  modulus) 
was  recorded  in  the  frequency  range  of  1  -100Hz.  These  tests  were  repeated  to  check  the  reproducibility 
of  the  results.  The  average  values  were  taken  for  both  the  runs  and  recorded  against  the  frequency  and 
are  presented  in  Figure  43  and  Figure  44. 
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Storage  modulus  vs  Frequency  (Hz) 
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Figure  43:  Storage  modulus  of  auxetic  inclusion  composites  tested  using  dual  cantilever  DMA. 
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Figure  44:  Tan  6  of  auxetic  inclusion  samples  acquired  using  dual  cantilever  DMA  measurements. 


These  results  shown  in  Figure  43  and  Figure  44  shows  a  generally  increasing  stiffness  with  increase  a- 
cristobalite  content  but  they  do  not  provide  a  clear  picture  of  the  influence  of  increasing  weight  fraction 
on  the  loss  behavior.  Specifically,  Figure  44  shows  that  the  largest  tan  6  values  were  observed  for  a  5% 
wt  composite,  while  the  lowest  values  were  observed  for  the  10%  wt  composite.  As  with  the  quasi-static 
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tests,  it  is  suspected  poor  bonding,  non-uniform  distribution  of  inclusions,  or  both  may  have  contributed 
to  this  observed  behavior.  Due  to  the  limited  time  of  this  project,  however,  a  conclusive  physical  basis 
for  the  observed  inconsistent  loss  behavior  was  not  able  to  be  identified. 

3.2.6.2  Modal  damping  analysis 

After  DMA  testing,  modal  analysis  tests  were  performed  at  the  University  of  Bristol  using  a  custom- 
made  rig  employing  electrodynamic  excitation  and  scanning  laser  vibrometry  (SLV).  The  tests  were 
carried  out  of  aluminium  beams  with  a  deposition  of  an  epoxy/a-cristobalite  composite  material 
patches  attached  to  the  beam  center,  as  shown  in  Figure  45a.  Typical  sample  dimensions  are  illustrated 
in  Figure  45a  and  the  test  setup  is  indicated  in  Figure  45b  and  c. 


Figure  45:  Details  of  the  force  vibration  measurements  of  modal  damping,  (a)  shows  the  dimensions  of  the  samples  (both  the 
aluminum  and  composite  patch),  (b)  shows  the  a  schematic  of  the  test  setup,  and  (c)  provides  a  photo  of  the  experimental 
setup. 


The  beams  were  initially  subjected  to  random  white  noise  to  identify  the  spacing  of  the  natural 
frequencies,  with  subsequent  pseudo-random  sweep  between  950  Hz  and  1050  Hz,  (corresponding  to 
the  3rd  flexural  mode  of  the  beam).  The  pseudo-random  sweep  was  driven  at  a  maximum  amplitude  of  1 
V  for  sweep  times  of  2  seconds.  The  input  force  was  measured  with  a  PCB  force  transducer  (calibration 
factor  42.24N/V,  with  range  of  200  mV).  Simultaneously,  a  SLV  (Polytec  PSV-400  F,  velocity  factor 
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25mm/s/V  with  channel  range  of  10V)  was  used  both  as  single  and  multi-point  velocity  acquisition 
device.  The  frequency  response  functions  (FRFs)  were  acquired  through  complex  averaging,  with  a 
central  frequency  of  1.02  kHz  and  0.1  kHz  of  bandwidth  and  201  spectral  lines. 

The  modal  damping  ratio  was  estimated  using  a  half-power  bandwidth  technique  implemented  in 
Modent  2008  software.  To  maximize  the  signal-to-noise  ratio,  the  technique  was  applied  to  the 
maximum  amplitude  point  of  the  FRF.  For  each  wt%  of  the  epoxy/amorphous  silica  coatings,  five 
samples  were  tested  to  acquire  a  sufficient  statistical  reproducibility  of  the  results.  Two  input  voltages  (1 
V  and  5  V)  were  considered  to  identify  any  nonlinearity  in  the  response.  The  mode  chosen  as  reference 
was  the  third  flexural  mode  because  it  resulted  in  the  highest  nodal  strain  for  the  free-free  configuration 
considered. 

An  initial  test  carried  out  on  a  reference  uncoated  aluminium  beam  has  provided  the  following 
values  for  the  undamped  natural  flexural  frequencies:  150  Hz,  500  Hz,  1000  Hz.  Figure  46  shows  the 
operational  mode  shapes  as  detected  by  the  SLV. 


First  Flexural  Mode 


Second  Flexural  Mode 


Third  Flexural  Mode 

Figure  46:  SLV  Images  of  the  first,  second,  and  third  flexural  modes  of  the  aluminum  beams  with  auxextic  inclusion  patches. 
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The  resulting  impact  on  the  patch  on  the  stiffness  of  the  beam-patch  composite  is  summarized  in  Figure 
47  which  shows  the  variation  of  the  third  natural  frequency  normalized  against  the  analogous  one 
belonging  to  the  uncoated  beam  sample. 
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Figure  47:  Variation  in  the  resonant  frequency  of  third  mode  of  vibration  for  two  different  driving  amplitudes  as  a  function  of 
the  weight  fraction  of  auxetic  inclusion  patches. 

These  experimental  results  show  that  there  is  no  apparent  dependence  of  the  natural  frequency 
behavior  against  the  value  of  the  input  force.  The  results  suggest  that  the  assumptions  of  linearity  of  the 
system  (at  least  in  terms  of  natural  frequencies)  are  satisfied.  However,  it  is  worth  noticing  the  large 
standard  deviations  existing  for  the  5  wt%,  20  wt%  and  25  wt%  samples.  It  is  apparent  that  a  more  than 
a  15%  increase  in  natural  frequency  occurs  for  the  15  wt%  samples.  The  assumption  of  linearity  is 
weaker,  however,  for  the  identification  of  the  modal  damping  ratio,  as  it  can  be  observed  in  Figure  7.  An 
average  20  %  increase  of  the  modal  damping  ratio  is  observed  passing  Frms  =  65  mN  to  Frms  =  923  mN. 
Due  to  the  amplitude  of  the  standard  deviations  associated  to  the  measurements,  it  is  difficult  to 
identify  a  clear  trend  existing  between  modal  damping  ratios  and  amount  of  amorphous  silica 
reinforcement.  However,  it  does  appear  that  a  clear  increase  in  loss  factor  exists  for  the  samples  with 
the  highest  wt%  considered  and  the  highest  amplitude  of  the  root  mean  square  force. 
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Figure  48:  Variation  in  the  modal  loss  factor  (in  percent)  of  third  mode  of  vibration  for  two  different  driving  amplitudes  as  a 
function  of  the  weight  fraction  of  auxetic  inclusion  patches. 

3.2.7.3  Ultrasonic  characterization 


The  final  set  of  characterization  testing  done  to  understand  the  behavior  of  the  auxetic  composite 
materials  was  ultrasonic  measurement  of  the  effective  stiffness  and  loss  in  a  frequency  band  around  1 
MHz.  The  test  was  accomplished  using  a  contact  ultrasonic  method  first  described  by  Trieber  and  co¬ 
workers  [44],  The  test  consists  of  six  separate  measurements  which  combine  reflected  wave  and 
transmitted  wave  acquisition  to  simultaneously  measure  the  phase  speed  and  attenuation.  These 
separate  measurements  permit  direct  measurement  of  the  reflection  coefficient  between  the 
transducer  and  the  material  sample  which,  unmeasured,  would  negatively  influence  the  accurate 
estimate  of  the  attenuation  coefficient  in  the  material.  The  parameters  directly  measured  during  the 
tests  are  the  phase  speed,  cph,  and  attenuation  coefficient,  a.  For  plane  wave  propagation,  these 
parameters  are  related  to  the  wave  disturbance  as  shown  in  Eq.  (30), 

u  ^z,  t )  =  uQe^ut  kz ^  =  uQe~aze^  ^ .  (30) 

In  this  formulation,  the  wavenumber  is  complex  and  defined  with  real  and  imaginary  parts  given  by 
k  =  k'  —  ja  .  The  attenuation  coefficient  is  a  metric  of  the  rate  at  which  a  forward  progressing  plane 
wave  amplitude  decreases  in  space  due  to  inherent  losses  in  the  medium.  The  real  part  of  the  wave 
number,  k’,  is  related  to  the  phase  speed,  cph,  through  the  relationship  k  =ujc ph  .  Using  the 

waveforms  captured  during  the  six  measurements  shown  in  Figure  49  it  is  possible  to  get  estimates  of 
the  stiffness  and  loss  properties  of  each  sample,  examples  of  which  are  shown  in  Figure  52  and  Figure 
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53.  During  ultrasonic  testing,  six  samples  of  each  weight  fraction  (0,  5,  10,  15,  20,  and  25%)  were 
measured  to  get  statistical  information  on  the  properties  of  the  composites. 
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Figure  49:  The  six  separate  measurements  associated  with  the  ultrasonic  characterization  method  used  to  characterize  the 
auxetic  inclusion  samples.  Measurements  Ml,  M2,  M5,  and  M6  are  use  reflected  signals,  while  measurements  M3  and  M4 
use  through  transmission  signals.  Combinations  of  those  signals  permit  direct  calculation  of  the  complex  valued  frequency 
dependent  reflection  coefficients  at  each  interface. 


The  tests  schematized  in  Figure  49  were  realized  by  creating  a  test  fixture  like  the  one  schematized  in 
Figure  50.  This  fixture  permits  simple  and  repeatable  attachment  of  the  contact  transducers  while 
minimally  influencing  the  position  of  the  sample  under  test. 


Clamping 

screw 


Rubber 


Figure  50:  Schematic  of  the  fixture  created  and  used  to  complete  ultrasonic  attenuation  measurements.  The  'fluid  layer'  is 
the  ultrasonic  complant  used  to  improve  the  coupling  between  the  transducer  and  auxetic  composite. 


The  ultrasonic  tests  were  accomplished  using  traditional  ultrasonic  measurement  components.  A 
Panametrics  square  wave  Pulse-Receiver  supplied  the  source  voltage  and  signal  conditioning  for  the 
receive  transducer  (whether  it  be  in  pulse-echo  or  transmit-receive  configurations).  The  transducers 
used  were  two  Panametrics  V603  Videoscan  1  MHz  center  frequency  longitudinal  reciprocal 
transducers,  the  signal  was  acquired  using  a  Tektronix  MS05000  series  digital  oscilloscope.  A 
representative  signal  acquired  for  Ml  (pulse-echo  with  a  free  surface)  is  shown  in  Figure  51.  As  was  the 
case  in  all  measurements,  a  post-processing  algorithm  extracted  the  first  and  second  arrivals  through 
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peak  detection  and  time  windowing.  The  resulting  spectrum  for  the  full,  first  and  second  arrivals  are 
shown  in  the  right  column  of  the  same  figure.  The  frequency  dependent  difference  in  magnitude 
between  the  first  and  second  arrival  is  indicative  of  the  signal  attenuation  (though  this  is  not  corrected 
for  reflection  losses  at  the  transducer-material  interface)  while  the  difference  in  phase  in  combination 
with  the  sample  thickness  can  be  used  to  deduce  the  phase  velocity  of  the  sample. 

Time  Domain  Signal,  Test 
Ml,  1st  Arrival 

5%  Sample  2  FFT  of  Received  Signal 


Time  Domain  Signal,  Test 
Ml,  2nd  Arrival 
5%  Sample  2 


Figure  51:  Representative  waveform  and  associated  spectrum  from  a  single  measurement  on  an  auxetic  inclusion  sample. 

As  an  example,  estimates  of  the  frequency  dependent  attenuation  coefficient  and  phase  speed  for 
various  samples  are  shown  in  Figure  52  and  Figure  53.  These  results  suggest  that  the  materials  show  a 
slight  anomalous  dispersion  (phase  speed  decreasing  with  increased  frequency)  and  increasing 
attenuation  with  increasing  frequency.  The  anomalous  dispersion  may  likely  results  from  processing 
artifacts,  especially  well  away  from  the  center  frequency  of  the  transducer,  1MHz.  The  increase  in 
attenuation  coefficient  with  increasing  frequency  is  typical  of  most  materials  an  due  to  a  combination  of 
material  losses  and  scattering  at  the  microscale. 
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Figure  52:  Measured  frequency  dependent  sound  speed  for  representative  material  samples  with  varying  weight  fraction  of 


a-cristobalite. 


Average  Attenuation  Coefficient, 
0%  Sample  2 


Average  Attenuation  Coefficient, 
5%  Sample  2 


Frequency,  (MHz) 


Final  Report 


New  Solutions  for  Energy  Absorbing  Materials 


Page  61  of  68 


Figure  53:  Attenuation  coefficient  measured  using  contact  ultrasonic  methods  of  representative  samples  for  varying  weight 


fractions  of  a-cristobalite  inclusions. 

Using  data  like  that  shown  in  Figure  52  and  Figure  53,  one  can  directly  the  measured  effective  stiffness 
and  loss  of  the  auxetic  composite  materials  using  Eqs.  (31)  and  (32): 


M  (/)  —  A  +  2/i 


1  +  Vl  +  772 

2(l +  772) 


(31) 


(32) 


The  results  of  these  calculations  for  all  samples  are  shown  in  Figure  54  which  plots  the  effective  plane 
wave  modulus  and  loss  factor  of  the  material  at  1  MFIz.  Note  that  the  raw  and  windowed  data  set 
shown  in  Figure  51  is  just  one  of  six  time-series  that  was  taken  for  a  single  sample  of  a  single  weight 
fraction  of  inclusion.  All  told,  the  entire  data  set  used  to  generate  estimates  of  the  stiffness  and  loss  of 
the  auxetic  composite,  shown  in  Figure  54,  consisted  of  36  measurements  per  sample,  or  216  individual 
measurements  to  measure  the  effective  frequency-dependent  properties  of  these  composite  materials. 
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Plane  Wave  Modulus  at  1  MHz 


Loss  Factor  at  1  MHz 


Figure  54:  Plane  wave  modulus  and  loss  factor  as  a  function  of  volume  fraction  a-critobalite.  Increases  in  both  stiffness 
(~30%)  and  loss  (~54%)  are  observed  for  the  auxetic  composite  10%  by  volume  composite  over  the  pure  matrix  material. 


Results  from  the  ultrasonic  test  seem  to  more  clearly  indicate  increases  in  stiffness  and  loss  as  a  function 
of  increasing  volume  fraction  of  auxetic  inclusions.  Notably,  the  stiffness  is  shown  to  increase  nearly  30% 
while  the  loss  factor  increases  almost  50%  with  the  use  of  only  10%  by  volume  of  auxetic  inclusion.  The 
augmentation  of  stiffness  for  this  amount  of  cristobalite  is  in  line  with  the  measurements  shown  in 
Figure  42,  while  the  enhancement  in  effective  loss  factor  is  greater  than  and  more  monotonic  than  what 
was  measured  using  either  DMA  or  vibration  techniques.  Overall,  the  ultrasonic  measurement  displays 
much  less  ambiguous  trends  than  the  results  of  the  DMA  and  vibration  measurements.  However,  due  to 
the  statistical  nature  of  these  results  and  the  clear  statistical  significance  of  these  trends,  we  are  certain 
that  the  ultrasonic  measurements  indicate  increases  in  lossy  behavior  of  auxetic  composites  at 
ultrasonic  frequencies. 

5  Insights  and  Future  Work 

The  results  presented  in  this  report  were  achieved  through  modeling  (analytical  and  finite  element) 
of  specific  structures  designed  to  exhibit  NS  and  experimental  characterization  of  NPR  inclusion 
composites.  The  work  associated  with  this  project  has  shown  that  structures  that  display  NS  behavior 
can  indeed  be  designed  using  buckled  elements.  Further,  a  very  general  finite  element  energy  based 
approach  has  been  derived,  implemented,  and  benchmarked  to  determine  the  nonlinear  effective 
stiffness  of  a  structured  element.  This  approach  is  unique  in  the  literature  and  very  powerful  for  the 
future  design  of  engineered  microstructure,  regardless  of  ultimate  application.  This  model  provides  a 
powerful  tool  in  the  design  of  microstructured  composite  materials.  This  project  also  integrated  this 
nonlinear  homogenization  model  with  more  traditional  multiscale  micromechanical  models  to  produce  a 
multiscale  model  that  predicts  macroscopic  stiffness  and  loss  parameters  as  a  function  of  microstructure 
and  mesoscopic  parameters  like  inclusion  fraction  and  orientation  distribution.  This  multiscale  model 
was  then  employed  to  do  an  exhaustive  design  space  exploration  of  a  family  of  structured  inclusion 
geometries.  That  exploration  provided  seed  data  to  generate  a  Kriging-based  surrogate  model  of  the 
influence  of  microscale  structure  on  mesoscale  anisotropic  stiffness.  This  combination  of  models 
provides  a  very  powerful  representation  of  this  nonlinear  multiscale  system  that  will  enable  rapid  design 
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of  these  materials  for  structure  damping  treatments  [45].  In  future  work,  NS  research  of  this  nature 
should  focus  on  the  implementation  of  the  models  validated  during  this  project  to  perform  an  iterative 
exploration  of  the  micro-  and  mesoscale  design  space  to  improve  the  elastic  and  absorptive  properties 
of  viscoelastic  composites.  Multiscale  design  model  should  be  wrapped  by  a  nonlinear  optimization 
scheme,  such  as  kernel-based  Bayesian  network  mappings  to  identify  optimal  microscale  designs  that 
lead  to  amplified  stiffness  and  loss  performance  through  the  use  of  NS  [46]. 

Auxetic  composite  material  work  associated  with  this  research  focused  on  the  fabrication  of  NPR 
composite  materials  and  their  stiffness  and  loss  performance.  It  provided  a  broadband  characterization 
of  the  viscoelastic  response  of  composites  containing  auxetic  inclusions,  ranging  from  static  to  ultrasonic 
frequencies.  This  was  achieved  through  tensile  tests,  dynamic  mechanical  analysis,  modal  damping 
analysis  with  damping  patches,  and  ultrasonic  measurements.  The  experimental  data  showed  mixed 
results  on  the  overall  stiffness  and  loss  of  the  composite.  Low  frequency  measurements  (DMA  and 
modal  analysis)  suggest  that  no  clear  correlation  between  increases  in  the  inherent  stiffness  and  loss 
properties  of  the  composite  with  increases  in  volume  fraction  of  auxetic  inclusions.  However,  ultrasonic 
measurements  unambiguously  indicated  that  both  stiffness  and  loss  increase  with  increases  in  auxetic 
fractions.  The  physical  root  causes  of  these  observations  were  not  clearly  determined.  The  latter  results 
suggest  that  auxetic  inclusion  composites  may  indeed  present  a  unique  means  to  engineer  high  loss 
composite  materials  without  degrading  the  overall  stiffness  of  the  material,  but  more  work  must  be 
done  to  validate  this  hypothesis.  Future  work  should  focus,  through  both  theoretical  and  experimental 
methods,  on  the  physical  phenomena  associated  with  energy  absorption  using  auxetic  inclusions. 
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